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^ ' Abstract 

In this paper, we mainly consider the problem of spherical distribution of 5 points, that is, how to 
configure 5 points on a sphere such that the mutual distance sum attains the maximum. It is conjectured 
that the sum of distances is maximal if 5 points form a bipyramid configuration in which case two points 
are positioned at two poles of the sphere and the other three are positioned uniformly on the equator. 
We study this problem using interval methods and related technics, and give a proof for the conjecture 
through computers in finite time. 
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1 Introduction 



Studies on the problem of optimally arranging points on a sphere can date back to over one hundred 
' years ago, when Thomson attempted to explain the periodic table in terms of the "plum pudding" model of 

fT^ , the atom. Since then, several varied problems were proposed, and some of such problems are still unsolved 

On ' now [1]. In general, these problems involve finding configurations of points on the surface of a sphere that 

maximize or minimize some given quantities, some of them are directly relevant to physics or chemistry 
' where stable configurations tend to minimize some form of energy expression. 

I The problem has the following general form. Let xi,X2, ■ ■ ■ ,Xn be points on the unit sphere 5"'^^ of the 

Euclidean space R™, denote 

y(X„,m,A)= (1-1) 

l<i<j<n 

' where Xn = {xi,X2, ■ ■ ■ ,a;„), and \xi — Xj \ denotes the Euclidean distance between Xi and xj. 



For A < 0, denote 



where 



Vi(n,TO,A)= min V{X,n,m,X) (1.2) 
Vi(n,m,0)== min V log r r (1.3) 

^,.CS"'-i ^ \Xt-Xj\ 

When TO = 3, this is the 7th Problem listed by Steve Smale in Mathematical Problems for the Next 
Century [2, •']]. 

For A > 0, denote 

V2(n,m,X)= max V(X,n,m,X) (1.4) 

So far as we know, G. Polya and G. Szego [4] first studied problems of such types in 1930s, since then, 
a number of results about V2(?t., to. A) have been derived. For example, L. Fejes Toth proved results for 
cases when to = 2, A = 1 and when ri, = TO + l,A = l[(i]. E. Hille considered the asymptotic properties of 
V2{n, TO, A)/A^ when n oo for definite to and A, and gave some results [7]. K. B. Stolarsky proved bounds 
of V2{n,m, X) for definite to and A in [8, 9], and gave some properties of point distributions corresponding 
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V2{n^ m, A) when m — 2 and m = 3 in [10, 11, 12]. R. Alexander also proved bounds of V2(?t., 3, 1) in [L!], and 
discussed some generalized sums of distances in [14, 15]. G. D. Chakerian and M. S. Klamkin proved bounds 
of l^(n,m, 1) in [1(5]. J. Berman and K. Hanes proved a property of the point distribution corresponding 
V2(n, 3, 1), and deduced some numerical results in [17]. G. Harman, J. Beck, T. Amdeberhan proved bounds 
of V2{n, m, A) in [18, 19, 24]. Similar problems were also discussed in [20, 21, 22, 23]. 

For V2(5,3, 1), numerical computations show evidences for the conjecture that, it is obtained when 5 
points form a bipyramid configuration in which case two points are on the two poles of S^, while three other 
points are uniformly distributed on the equator. In this paper, we study this problem via interval arithmetic, 
and prove the conjecture through computer in comparatively short time. For related problems, this guides 
a different method. 

The main ideal of our proof is as follows. Firstly we express V{X^,3^1) as a function under certain 
coordinate system, secondly we exclude a domain where the bipyramid configuration is proved to correspond 
an only maximum of ^(Xs, 3,1), lastly we subdivide the remaining domain, and prove that function values in 
these subdomains are less then the previous maximum obtained. So we complete the proof of the conjecture. 

2 Mathematical descriptions of the problem 

2.1 Spherical coordinate system 

We choose the spherical coordinate system as showed in Fig. 1. A point P on is identified by (1, cj), 9), 
where G [~f ' f ] ^^'•^ angle from vector OH, i.e., the projection of vector OP in xoy-plane, to vector 
OP, positive if the z-coordinatc of P is positive, and 9 e [— tt, tt) is the angle from x-axis to vector OH, 
positive if the y-coordinate of P is positive. 



According to such definitions, we have following formulas transforming from spherical coordinate (1, (j), 9) 
to Cartesian coordinate (x, y, z), 




y 



Fig. 1: The spherical coordinate system 




cos(^) cos(0), 
cos(^) sin(6'), 
sin ((/)). 



(2.1) 
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Considering the spherical symmetry, we can choose the spherical coordinates for 5 points as 



follows: 



^(1, 0, 0), 5(1, 01, tt), 02, ^2), D{1, 03 


03),S(l,04,e 


4). 


sum of mutual distances of these points is 












— y ^ + z cos [(pi) + \j I — I cos ((p2 ) cos [pi ) 






+ V2 - 2 cos (03) cos (6*3) + a/2 - 2 cos (<, 


64) COS (04 ) 




+ \/2 cos (0i) cos (02) COS (^2) + 2 — 2 sin 


(01 ) sin (02) 




+ \/2 COS (0i) cos (03) cos (^3) + 2 — 2 sin 


(0i)sin(03) 




+ -\/2 cos (0i) cos (04) cos (^4) + 2 — 2 sin 


(0i)sin(04) 




+ cos (03) COS (02) COS (6*2 - 6I3) + 2 - 


- 2 sin (02) sin 




+ V-2 COS (02) COS (04) COS (6*2 - 6I4) + 2 - 


- 2 sin (02) sin 


('/>4) 


+ V — 2 COS (03) COS (04) COS (6*3 — 64) + 2 - 


- 2 sin (03) sin 


(</>4) 



(2.2) 



(2.3) 



2.2 Bipyramid distribution 

Spherical coordinates of 5 points corresponding a bipyramid distribution are not unique, but the following 
5 points indeed form a bipyramid configuration. 



A(l,0,0),i?(l,-|,7r),C(l,|,^),D(l,0,-|),£;(l,0,|), 



as showed in Fig. 2. 



(2.4) 




Fig. 2: The bipyramid distribution 
Denote the corresponding values of (0i, 02, 03, ^'a, 04, 6*4) by 
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then the corresponding value of function / is 

fmax 

and the Hessian matrix of / is 



3\/3 + 6V2 + 2 
15.68143380, 



(2.5) 
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(2.6) 
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This matrix is negative definite, so the bipyramid distribution corresponds a maximum of function /. 

2.3 Inequality form 

As a matter of fact, what we are to prove is the fohowing inequahty, 

02, 6*2, 03, 6*3, 04, Oi) < fmax, (0i, 02, 6*2, 03, 6*3, 04, 0i) e v. 



(2.7) 



-7r,7r), 



where 

V^"2'2^'' 2'2''' 2'2''' 2'2'" 

and the equahty holds if and only if (0i, 02, ^2, 03, ^3, 04, ^4) = ©fcp- 

In the remaining part of this paper, we will according to following steps to prove this inequality. 



TT TT TT TT 

--,-], h^,7r),[--,-] 



1. Giving some restricted conditions and results to demonstrate that we only need to prove the inequality 
over a subdomain of V, i.e., P^^' U P^^) (ggg gq. (4.6)). 

2. Analyzing interval Hessian matrices (Theorem 4.3, 4.4 and 4.6) to prove that the equality holds only 
at Qbp over a subdomain of V^^^ U V^^^ i.e., V^p (see Proposition 4.1). 

3. Analyzing interval Hessian matrices (Theorem 4.5 and 4.7) to prove the corresponding strict inequality 
holds over a subdomain of V^^^ U V^^^ i.e., Vp (see Proposition 5.1). 

4. Making use of the interval arithmetic (§ 4.1.1) to prove the corresponding strict inequality holds over 
the remaining domains, i.e., (V^^^ UV^'^^)\{VbpUT>p) (see Eq. (6.1)). 



3 Restricted conditions and verification domain 
3.1 Some results 

What we are to prove is in fact that, there exists no distribution of 5 points exclude the bipyramid 
distribution corresponding larger distance sum then fmax. We need following results so as to simplify this 
problem. 
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Proposition 3.1. If some configuration of 5 points corresponds larger function value of f then the bipyramid 
configuration, and AB is the second largest distance in (2) ~ 10 distances, 0i should satisfy 

(j)i > -2 arccos(\/3/6 + V2/3) (3.1) 

Proof. From Equation (2.5), we know that in order to attain larger distance sum than that the bipyramid 
configuration corresponds, the second largest distance must be not smaller then 

((3 V3 + 6 V2 + 2) - 2)/9 = \/3/3 + 2 V2/3. 

With the condition that AB is the second largest distance, the result required can be deduced immediately. 

□ 

Proposition 3.2. // 5 points are on the same half sphere, f can not attains its maximum. 

Proof. Without loss of generality, suppose z-coordinates of 5 points are all nonpositive, if the 2-coordinate 
of some point is negative, we move it to the symmetric position with respect to the xoy-plane, then we will 
get a larger distance sum. 

If 5 points are all distributed on the xoy-plane, the maximal distance sum is [li] (5 points form a regular 
pentagon) 5 cot , which is obviously smaller then the mutual distance sum corresponding the bipyramid 
configuration (see § 2.2). □ 

Proposition 3.3. // a partial derivative of function f does not vary signs in a domain, then there exists no 
stationary point of f in this domain. 

Theorem 3.1. [17] Let pi, . . . ,p„ be points on the unit sphere S'^ in R^. Let f : S'^ ^ M. be defined by 

n n 

/(^) = X] 1^ ~ -Pil ■ ^/ / ^^'5 maximum at p, then p = (j/l^l, where q ^ (-P ~ Pi) / \P ^ Pi\- 

4=1 4=1 

Theorem 3.2. [S] Suppose the 5 points are placed so that function f is maximal, then any distance between 
two points cannot be less then j^. 

3.2 Some restricted conditions 

We can consider the problem under following restricted conditions due to above results. 
Condition 3.1. AB is the second largest distance in (2) = 10 distances. 

Condition 3.2. D is on the left half sphere, C,E are on the right half sphere, C is above E. 
Condition 3.3 (by Proposition 3.1). (jyi > —2 arccos(V3/6 + a/2/3) 
Condition 3.4 (by Proposition 3.2). Five points are not on any half sphere. 
Condition 3.5 (by Theorem 3.2). Distances between any two points are larger then j^. 

3.3 Domain subdivision 

Under these conditions, the bipyramid configuration (corresponding the maximal distance sum con- 
jectured) and the pyramid configuration (corresponding another stationary point of the function /) each 
corresponds only one coordinate representation. Further more, we can divide the domain in which we need 
to verify no distribution of points corresponds a larger distance sum into the following two subdomains: 

1. D is on the upper half sphere (denote this domain by V^^^): 

01 G [-2 arccos (^/3/6 + V2/3) ,0], 

02 e [-7r/2,0], 

03 € [0,7r/2], 

03 e h7r,0], 

04 € [-7r/2,0], 
04 e [0,7r]. 
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2. D is on the lower half sphere, C is on the upper half sphere (denote this domain by P'^^): 

01 e [-2 arccos (V3/6 + V2/3) ,0], 

02 e [0,7r/2], 
^2e [0,7r], 

03 e h7r/2,^/2], 
Ose [-7r,0], 

04 e h7r/2,7r/2], 
04 G [0,7r]. 

Now, we are to prove that, under Condition 3.1 - 3.5, function / attains its maximum in and 2?(2) 
at the only point corresponding the bipyramid distribution of A, B, C, D, E, i.e., 

/(01, 02, 02, 03, 03, 04, 04) < fmax, (01, 02, 02, 03, 03, 04, 04) G V'^'^ U P^^). (3.2) 

where the equality holds if and only if (0i, 02, 02, 03, 03, 04, 04) = Qbp- 

In the following parts of this paper, we will illustrate the domain verification methods, and detailed steps 
as well as results. 

4 Domain near coordinates corresponding the bipyramid distri- 
bution 

4.1 Interval methods 

We first briefly introduce the interval methods we used in our proof. 

4.1.1 Interval arithmetic 

We define an interval as a set [25]: 

X :^[a,b]^{x:a<x<b}, (4.1) 

where a, 6 e M. 2^,X respectively denote the left and right vertexes of the interval X. 

For intervals X and Y, ii x > y for each x € X and each y G Y, we say that X > Y. Other interval 
relations are understood the same way. An 7i-dimensional "interval vector" is an n-tuple of intervals X = 
(Xi, . . . ,Xn), which is used to denote some rectangular domain in R". Let IM be the set of intervals over 
K, and be the set of n-dimensional "interval vectors" . 

We can define an imbedding from M to IM as follows 

/i(a;) = [x,x], 

thus for numbers in R, we can also consider them as intervals. 
We define interval arithmetic over IR as 

XoY^{xoy: xeX,ye Y}, 

where o is " + "," — ","* " or "/". Further more, for an elementary function /, we define a corresponding 
elementary mapping as 

fix) = {fix) -.xex}. 

When operands of interval arithmetic or arguments of elementary functions are intervals, we consider 
underlying computations are interval computations defined above, and the interval computation is of the 
same precedence as the corresponding arithmetic computation. 

Under above definitions, an arbitrary elementary function / ; R" R can be expanded to a mapping 
over IR" IR: 

/(X) = /(X). (4.2) 
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Through such /, We can get an interval which contains the function range of / over rectangular domain X, 
this is the critical point we solve the problem. 

As a matter of fact, there are related programs used to process interval arithmetic, such as the procedure 
evalr can be used to implement interval arithmetic without errors. But in practice, it may be not necessary 
to implement errorless interval arithmetic, because what we get from interval arithmetic are just intervals 
contains the ranges of function values. Another problem is that acting such errorless interval arithmetic is 
always time-consuming, thus it cannot meet our needs. 

Considering the efficiency and the accuracy, we wrote an interval arithmetic package IntervalArith- 
metic based on the Maple system. The package uses rational numbers as interval vertexes, and acts com- 
putations with controllable errors. In fact the result it computes for /(X) is an larger interval containing 
/(X), and the difference can reduce to zero as intervals of X shrink to points. For the detailed code, see 
Appendix A.l. 



4.1.2 Interval matrices 

Relations of real matrices of the same order are understood componentwise. An interval matrix is defined 
as the following set of matrices: 

([%, ^[A,A] = {AeW^ ■.A<A< A}, 

where 

A^{a,^),A^ia,j). 

When A and A are symmetric, we call the set of symmetric matrices in [A, A] a symmetric interval matrix 

which is also denoted by [A, A] . 

— A + 'A 'A — A 

For a interval matrix [A, A] , denote its midpoint matrix by Ac — — , radius matrix by Ag — — . 

For a real symmetric matrix A, it is well know that all its eigenvalues are real, we denote them in decreasing 
order by Xi{A) > X2{A) > ••• > A„(A), and denote the spectrum of A (i.e. the maximum eigenvalue 
modulus) by g{A). For bounds of eigenvalues of matrices in an interval matrix, it can be directly deduced 
from the Wielandt-Hoffman theorem [28] that 

Theorem 4.1. For a symmetric interval matrix [A^A], the set 

{XM)--Ae[AA]} 
is a compact interval, denote this compact interval by 

[X,{[A,A])M[A,A])],'^<i<n, 

then 

M[A,A])M[AA])] C [X^iAc) - g{As), X,{Ac) + giAs)], i^l,...,n. 
In fact. Ai([A, y4]) and A„([A, A]) can be solved explicitly [2()], that is, 
Theorem 4.2. A real symmetric interval matrix 

([a.,, a,,]) = {A e M" " : A < A < = (a,,), A = (ay)} 
corresponds following 2"^^ vertex matrices: 

Afe = (afe,,),0<fc<2"-i-l, 
where we denote the binary representation for k by k = (fcifc2 ■ ■ ■ fc„)2, and 

= ^K- +a,j + (-l)'=-+^^(ay- - ay)). 

For matrices in this symmetric interval matrix, minimal ( or maximal) eigenvalues of them attain the mini- 
mum (respectively, maximum) at some vertex matrix A^. 
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For a real symmetric interval matrix [A, A] , we say it is positive (semi)definite if A is positive (semi)definite 
for each A G [A, A], and it is nonpositive (semi)definite if A is not positive (semi)definite for each A E [A, A]. 
Definitions such as negative (semi)definiteness, nonnegative (semi)definiteness of [A, A] are understood the 
similar way. 

Now we introduce the results for verifying positive definiteness and nonpositive semidefiniteness of sym- 
metric interval matrices, which can directly deduce criterions for negative definiteness, nonnegative definite- 
ness, etc. 

Rohn has given the following theorem [26], which is an improvement on results in [27], we state it in a 
varied way that adapts to be understood as an algorithm. 

Theorem 4.3. The real symmetric interval matrix 

([a.,, ay]) = {A G M" X " : A < A < = {a.,^),A= (Sy)} 
is positive definite if and only if the following 2"~^ vertex matrices are all positive definite: 

Afc = (afe,,),0<fc<2"-i-l, 
where we denote the binary representation for k by k ^ (^1^2 ■ • • kn)2, and 

aMj = +a,j + (-l)'=-+^^(ay - ay)). 

Theorem 4.2 in fact implies Theorem 4.3, because the positive definiteness of a symmetric interval matrix 
is equal that the minimum of minimal eigenvalues of matrices in it is positive. So we can get from Theorem 
4.1 a sufficient condition for determining the positive definiteness of a symmetric interval matrix. 

Theorem 4.4. The symmetric interval matrix [A, A] is positive definite if 

Xn{Ac) - ei^s) > 0, 
where Ac and As are the midpoint matrix and the radius matrix respectively. 

Similarly, the nonpositive semidefiniteness of a symmetric interval matrix is equal that the maximum of 
minimal eigenvalues of matrices in it is negative, i.e.. 

Theorem 4.5. The symmetric interval matrix [A, A] is nonpositive semidefinite if 

Xn{Ac) + Q{As) <0, 

where Ac and Ag are the midpoint matrix and the radius matrix respectively. 

The procedure isdef in the package IntervalArithmetic can implement algorithms above to verify the 
properties of symmetric interval matrices, such as positive definiteness, negative definiteness and so on. 

With the help of above theorems, we can use the following trivial results to determine the extreme point 
of a function in a domain. 

Theorem 4.6. If g G C'^{D), Xo G D is a stationary point of g, and the Hessian matrix of g over D varies 
in a positive definite real symmetric interval matrix, then Xq is the minimum point of g in D. 

Theorem 4.7. If g G C^{D), and the Hessian matrix of g over D varies in a nonpositive semidefinite real 
symmetric interval matrix, then no inner point of D is the minimum point of g. 
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4.2 Domain excluded near coordinates corresponding the bipyramid distribu- 
tion 

Now we introduce a disturbance [— gfy, gfy] on coordinates corresponding the bipyramid distribution, 
and obtain a rectangular domain, i.e., 



In this domain, 



(^2 

03 



V ^4 ) 



( [--380- TT --2Ii-7r] \ 

/ L 1131 1131 "J \ 



131 ' 1131 J 



[376 378 , 



L377 



377 



377 



377 

-379 ^ - 375 1 
754 " ' 754 "J 



V 



377 

"375 
.754 ^ 



' 377 

379 „ 
754 ^ 



(4.3) 



92 varies in tt, |^ tt] , which exceeds the bound we prescribed for 6*2. But since the 
periodicity of function /, it is of no error. In fact, interval vertexes are represented by rational numbers in 
the Maple package IntervalArithmetic, so these intervals whose vertexes contain tt are enlarged to their 
rational representations, that is, the rectangular domain we actually obtain is 



hp 



1055530689 
' 1000000000 ' 



1038864413 1 
1000000000 J 



[ 1038864413 1055530689 1 
L 1000000000' 1000000000 J 

[ 783314879 98435181 1 
L 250000000 ' 31250000 J 



10416421265218147343 



10416421265218147343 



1250000000000000000000 ' 1250000000000000000000 J 



315825893 
200000000 ' 



10416421265218147343 



1562463189 1 
1000000000 J 



10416421265218147343 



V 



1250000000000000000000 ' 1250000000000000000000 J 

[1562463189 315825893 1 
L 1000000000 ' 200000000 J 



(4.4) 



The interval Hessian matrix of / over TDhp can be calculated by interval arithmetic: 

V = (1/1,^2, V^3,"^^4,"t1s,V^6,VV), 

where Vi(z = 1, . . . , 7) are vectors as follows: 



(4.5) 



r 9073071021 
L 10000000000 
r 4158136493 
L 10000000000 ' 

2503191333 
1000000000000 ' 
r 910889963 
L 2500000000 ' 
r 6038013477 
L 10000000000 ' 
r 72871197 



4158136493 1126402921] 



200000000 
3104766243 



312500000 J 
1126402921 -I 
2500000000 J 

2503191333 l 
1000000000000 J 
428510269 1 
1250000000 J 
1552383121 1 
2500000000 J 
685616431 l 
2000000000 J 
6038013477 1 



L 10000000000 
2283824451 
2500000000 ' 
23364423 
- 781250000 ' 
114837637 
312500000 ' 
1559063347 
2500000000 ' 
183740219 
500000000 
r 6011315931 



2500000000 J 
8191611957 -I 
10000000000 J 
747661531 1 
25000000000 J 
3397366803 1 
10000000000 J 
240452637-1 
400000000 J 
424670851 1 
1250000000 J 
779531673 I 



5000000000 10000000000 J 



L 10000000000 1250000000 J 
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2503191333 
1000000000000 ' 
r 23364423 
L 781250000 ' 
3523502821 
10000000000 ' 
r 1628135751 
L 10000000000 ' 
r 1556247929 
L 20000000000 ' 
r 718117527 
L 5000000000 ' 
r 7781239691 
L 100000000000 ' 25000000000 J 

r 6038013477 15523831211 
L 10000000000 ' 2500000000 J 
r 1559063347 2404526371 
L 2500000000' 400000000 J 
r 1556247929 9916175537 l 
L 20000000000 ' 100000000000 J 
r 483991657 97163283 l 
L 20000000000 ' 4000000000 J 
r 1058713931 1002296947-] 



2503191333 1 
1000000000000 J 
747661531 1 
25000000000 J 
2901860369 1 
10000000000 J 
359058763 l 
2500000000 J 
9916175537 1 
100000000000 J 
1628135749 1 
10000000000 J 
2479043873 1 



. V4 = 



1000000000 
1041678267 



1000000000 J 
1041678267 



10000000000000 10000000000000 J 
[■312434903 12501736191 
1-626000000' 2500000000 J 

r 3104766243 6038013477 I 



L 5000000000 
r 6011315931 
L 10000000000 ' 
r 7781239691 
L 100000000000 ' 

1041678267 
10000000000000 
r 312434903 
L 626000000 ' 
r 607208011 
L 26000000000 ' 
r 1068713929 
L 1000000000 ' 



10000000000 J 
779531673 1 
1250000000 J 
2479043873 1 
25000000000 J 

1041678267 1 
10000000000000 J 
1250173619-1 
2500000000 J 
302494783 -| 
12500000000 J 
20046939-1 



2600000000 
114837637 
312600000 ' 
1628136761 
10000000000 ' 
214384921 
200000000 ' 
483991657 
- 20000000000 
r 10001389 
L 20000000 ' 
1041678267 



1250000000 J 
3397366803 -| 
10000000000 J 
359058763 l 
2500000000 J 
9891691177 l 
10000000000 J 
97153283 1 
' 4000000000 J 
4999306693 -| 
10000000000 J 
1041678267 



10000000000000 10000000000000 J 



L 200000000 
r 183740219 
L 500000000 ' 
r 718117527 
I- 6000000000 ' 
r 10001389 
1. 20000000 ' 
1041678267 
10000000000000 ' 
r 1071924603 
L 1000000000 ' 
r 607208011 



2000000000 J 
424670851 1 
1250000000 J 
1628136749 1 
10000000000 J 
4999305593 l 
10000000000 J 

1041678267 l 
10000000000000 J 
24729227991 
2500000000 J 
302494783 1 



25000000000 12500000000 J 



Through Theorem 4.3, we can judge that the symmetric interval matrix V is negative definite, and by 
Theorem 4.6, the conjectured configuration indeed corresponds the maximum of / in Vbp- That is 

Proposition 4.1. The bipyramid distribution of 5 points represented by Eq. (2.4) is the only maximal 
distance sum distribution in domain 'Dbp, 



/((^l, (^2, 6*2, 03, 6*3, 04, ^'4) < fmax, {4>1, (1)2,02,(1)3,63,41^,04) e Vbp, 

where the equality holds if and only if [(pi, (f)2, O2, 03, ^3, 04, O4) = Qbp- 



(4.6) 



5 Domain near coordinates corresponding the pyramid distribu- 
tion 

Under conditions in § 3.2, coordinates representing the pyramid distribution are unique, while they corre- 
sponds a stationary point of function /, and the function value on this point is too close to fmax, therefore, 
we discuss it separately. 

5.1 Pyramid distribution 

The spherical coordinate corresponding the pyramid distribution is 
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where 



^(1,0,0), 
5(1,-2^^1,^), 

D (1,^2,-^^3) , 
E{1, 0:2,(^3) , 



3 \/41 - 28 V2 

ixii = arcsin I 1 1 

'42 4 



a;2 = — arcsin 



3 V2 ^41 - 28^2^ 

4 ~2~ 



/ 3 V2 V4I " 28 ^2 V ^ 



3 V2 ^41 - 28 V2 
'4 ~2~ 4 



a;3 = arccot 



as showed in Fig. 3. 



4 2 



3 ^41 - 28 



4 2 



2\ 



(5.1) 




Fig. 3: The pyramid configuration 
Denote the corresponding vahies of (0i, 02, ^2, 037 ^3, 04, ^4) by 

Op = (^-2^1,^ -Wi,7r, 0^2,-^3,^2,^3^ , 

then the corresponding value of function / is 

pyfmax = /(6p) 



(5.2) 

15.67482117. ^ ^ 
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5.2 Domain excluded near coordinates corresponding the pyramid distribution 

Similarly with the method we adopted near the bipyramid distribution, we introduce a disturbance of 
[— on coordinates corresponding the pyramid distribution, and finally obtain a rectangular domain 



V 



5157880419 
10000000000 ' 



203137879 ] 
400000000 J 



' 1310916469 2637719631 



.1000000000' 200000000 J 
'156881049 31455643271 



50000000 ' 1000000000 J 



39276321 
156250000' 



1508635943 
1000000000 ' 



39276321 
156250000' 



2434251099 1 
10000000000 J 



375173149 1 
250000000 J 



2434251099 1 
10000000000 J 



L 250000000' 1000000000 J 



\ 



(5.3) 



The interval Hessian matrix of / over Vp can be calculated by interval arithmetic, i.e. 



(5.4) 



where Wi{i = 1, . . . ,7) are vectors 

/ r 68555671 80S5653SS7 "I , 



Wi = 



SOOOOOOO 
[241020403 



10000000000 J 
4060471777 1 



L 625000000 10000000000 J 
5396514777 5396514777 n 

10000000000000 ' 10000000000000 J 



1000000000 
r 323160259 
L 1250000000 ' 
6654530721 
10000000000 ' 
r 545472247 
1. 2000000000 ' 



3261193771 "I 
5000000000 J 
2727361243 1 
10000000000 J 
6522387553 l 
10000000000 J 
161580131 



, W2 = 



[241020403 4060471777 " 



L 625000000 
r 567279649 
L 500000000 ' 
609173849 
50000000000 ' 
68583739 
400000000 ' 
1486861643 
2500000000 ' 
428648367 
2500000000 ' 
[2937974709 



10000000000 J 
544204613 1 
500000000 J 
1218347717 l 
100000000000 J 
1584190467 1 
10000000000 J 
5875949419 1 
10000000000 J 
792095233 l 
5000000000 J 
5947446571 l 



L 5000000000 10000000000 J 



5396514777 



5396514777 



665453073 3261193771-] 



10000000000000 10000000000000 J 
r 609173849 1218347717 i 
L 50000000000' 100000000000 J 
r S0S5S54539 1541389841 -| 



100000000000 
9959577971 
100000000000 ' 
r 116759761 



25000000000 J 
9385595363 l 
100000000000 J 
109700733 1 



L 5000000000 4000000000 J 

- 0385595361 9959577971 1 

- 100000000000 ' 100000000000 J 

- 2335195223 2742518293 1 

- 100000000000 ' 100000000000 J 

r 323160259 2727361243 1 
L 1250000000 ' 10000000000 J 
1486861643 5875949419 1 



2500000000 
r 116759761 



10000000000 J 
109700733 1 



L 5000000000 4000000000 J 
585736263 2143830967 



- 5000000000 
r 36387877 
L 31250000 ' 
1002395789 
100000000000 ' 
r 4813603601 



25000000000 J 
225679327-1 
200000000 J 

616963279 1 
100000000000 J 
1215175283-1 



, Wo = 



L 1000000000 
r 68583739 
I. 400000000 ' 

9959577971 
100000000000 ' 

- 8S2S39S133 

- 10000000000 ' 
585736263 

L 5000000000 ' 
4897161727 
10000000000 ' 

r 616963279 



5000000000 J 
1584190467 T 
10000000000 J 
9385595363 l 
100000000000 J 
1678302277-1 
2000000000 J 
2143830967 -| 
25000000000 J 
4822526949 l 
10000000000 J 
1002395789 1 



L 100000000000 100000000000 J 



6654530721 
10000000000 ' 
r 428648367 
L 2500000000 ' 

- 9385595361 

- 100000000000 ' 

4897161727 
10000000000 ' 
1002395789 
100000000000 ' 
SS2S39S113 
. 10000000000 
r 171506479 



6522387553 1 
10000000000 J 
792095233 1 
5000000000 J 
9959577971 l 
100000000000 J 
4822526949 1 
10000000000 J 

616963279 1 
100000000000 J 
4195755701 l 
5000000000 J 
1171472521 -I 



L 10000000000 2500000000 J 



L 2000000000 10000000000 J 
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L 2000000000 62500000 J 
r 2937974709 5947446571 "I 
L 5000000000 ' 10000000000 J 

- 2335195223 2742518293 1 

- 100000000000 ' 100000000000 J 

- 616963279 1002395789 1 
. 100000000000 ' 100000000000 J 

r 4813603601 12151752831 
L 10000000000 ' 2500000000 J 
r 171506479 1171472521 n 
L 2000000000' 10000000000 J 
r 582206031 5641983191 



Through Theorem 4.5, we can judge that W is nonnegative semidefinite, and when the disturbance 
enlarges very httle, it is stiU true. So by Theorem 4.7, we know that values of / cannot attain the maximum 
in Vp, i.e.. 

Proposition 5.1. Maximum of function f can not attain in domain Dp, i.e., 

f{(j>i,(j>2,92,(l>3,03,(j)4,04)<fmax, (01, 02, 6*2, 03, 6*3, <A4, 6*4) e Pp. (5.5) 

6 Other domains 

Now, we are to prove the following strict inequality, 



/(0i, 02, 6*2, 03,^*3, 04, 6*4) < fmax, 

where (0i, 02, ^2, 03, ^3, 04, ^4) e (P'^' U V'-^^)\(Vbp U Vp). 



(6.1) 



Algorithms in this section are implemented by procedures in the Maple package fivepoints, for the code, 
see appendix. 

6.1 Branch and bound strategies 

We check domains over which variables take using the interval method, more precisely, we compute the 
interval value of the interval mapping corresponding some functions through interval arithmetic, properties 
of this interval may suggest that, when variables take values in this domain, function / has no stationary 
point, or its maximum is less than the value corresponding the bipyramid configuration, or it is not necessary 
to consider the case for symmetries. All in all, function values in this domain cannot be greater than fmax 
(these verification methods are implemented by procedure ischecked in the package fivepoints). The 
foUowings arc methods we used to exclude domains contained in V'-^^ and V^^y 

1. (by Condition 3.2) Verify that C is below E. 

2. (by Condition 3.4) Verify that 5 points arc in the same half sphere. 

3. (by Condition 3.1) Verify that AB is not the second largest distance. 

2 

4. (by Condition 3.5) Verify that the distance between some two points is less than — . 

15 

5. Verify that the upper bound of function values are less than fmax (see Eq. (2.5)). 

6. (by Proposition 3.3) Verify that some partial derivative of / does not change signs in this domain. 

7. Compute the interval Hessian matrix corresponding this domain, and determine its negative dcfinitcness 
through Theorem 4.3. 

8. Compute the interval Hessian matrix corresponding this domain, and determine its negative definiteness 
through Theorem 4.4. 
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9. Compute the interval Hessian matrix corresponding this domain, and determine its nonnegative defi- 
niteness through Theorem 4.5. 

10. Determine there exists no maximal point in this domain through Theorem 3.1. 

Different methods should be used in different domains, for example, methods 7 and 8 should be used 
first near points corresponding bipyramid distribution, then others (10, 5, 6) can be used; while method 9 
should be used first near points corresponding pyramid distribution, then others (10, 5, 6) can be used; and 
for generic domains, methods can be used in turns as 1, 2, 4, 3, 5, 10, 6. 

For a domain to be verified, we choose appropriate verification methods and the verification order, if 
verifications are not successful, we subdivide the interval whose width is maximal into two equal intervals, 
and verify the two subdomains recursively. We set a positive number, if the largest interval width of a 
domain we get in the above process is less than this number, we stop subdividing this domain, and record 
it, this domain may contain distributions of points corresponding larger distance sums then the maximal 
distance sum conjectured. This process terminates when all domains have been verified. If all domains are 
verified successful, and no domain is contained in the record list, then we have proved the conjecture in 
fact. The complete algorithm is described below (implemented as the procedure spchecked in the package 
flvepoints): 

Algorithm 1: CheckDomain 

Input: intcrvallist, methods, notcheckbipyiamid, notchcckpyramid 
Output: true/false 

1 begin 

2 chcckbipyramid :— not notchcckbipyramid; chcckpyramid :— not notchcckpyramid; chcckmethods :— methods 

3 if chcckbipyramid then 

4 if intcrvallist is contained in domain bipyramidintcrvallist then 

5 ^bipyramidintcrvallist is the rectangular domain we excluded first near the bipyramid distribution. 

6 add [—1] to checkprocess 

7 ^record the checking process 

8 return true 



9 
10 

11 
12 
13 
14 

15 
16 

17 
18 
19 

20 
21 
22 
23 
24 

25 
26 
27 
28 
29 
30 
31 

32 
33 

34 
35 



36 
37 end 



if some variable interval in intcrvallist disjoints the corresponding variable interval in bipyramidintcrvallist then 
chcckbipyramid :— false 

if chcckpyramid then 

if intcrvallist is contained in domain pyramidintervallist then 

i^pyramidintervallist is the rectangular domain we excluded first near the pyramid distribution. 
add [0] to checkprocess; return true 

if some variable interval in intcrvallist disjoints the corresponding variable interval in pyramidintervallist then 
chcckpyramid :— false 

while methodsi G methods do 

if checking domain intcrvallist successfully through the method methodsi then 
add [i] to checkprocess; return true 

:— the widest interval position in interval vector intcrvallist 
if the width of the dim-th interval in intcrvallist < ^qqq then 
add intcrvallist to notchecked 

^notchecked records domains cannot be checked successfully 
add [—4] to checkprocess; return false 

add dim to checkprocess 

^dim designates the variable to be subdivided, and is to be recorded in checking process 
subintcrvallist :— subdivide the domain intcrvallist in the dim-th interval 
cur :— true 

forall subintervallisti ^ subintcrvallist do 

if not CheckDomain(subintervallisti , chcckmethods, not chcckbipyramid, not chcckpyramid) then 
cur :— false 

if cur then 

add [—2] to checkprocess 

else 

add [—3] to checkprocess 
return cur 



6.2 Verification process 

In order to subdivide domains into appropriate widths, we act some experiments first, finally we subdivide 
domains as follows: for I?'^^ and I?'^' we divide in § 2, each domain is subdivided the way that each interval 
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of it is trisected, so we get 3^ — 2817 subdomains each, denote them respectively by: 



and 



7-,(l) 7-,(l) 7-,(l) 
'-^\ 1 '-^2 7 ■ ■ ■ I ■'^2187' 



If some of these subdomains are difficult to verify successfully, we can again subdivide them the same way. 
Actually, the following domains need to subdivide again: 

7-,(l) -nil) 7-,(l) -0(1) 7-,(l) 7-,{l) -rp) 

■^^62 ' ■'^ISS' ■'^239' ■'^SeS' ■'^1102! ■'^IIOS' '^1106' '^2114' ■'^2132' 

7-,(l) 7-,(l) 
-^1105-1101' -^1106-834' -^1106-861' -^1106-1099' -^llOe-llOO' 

7^(1) 7^(1) ^ ■ 

-^1105-1101-1100' -^1106-834-725' -^1106-834-726' 

7^(1) 

-^^1106-834-725-1752' -^^1106-834-726-1507' -^^1106-834-726-1750' 

where 2?J^Qg_j^]^Q]^ denotes the 1101-th subdomain in all 2187 subdomains of I?iio5, other similar notations 
are understood the same way. 

6.3 Algorithm implementations 

The Maple Package fivepoints implements algorithms described in above sections. For the detailed 
code, see Appendix A. 2. 

7 Conclusion 

The following is verification time for various domains (may differ on different computers, it is the time 
used by computers with Pentium IV 3.0 GHz CPU, and 1 GB RAM): 

1. Time used to verify domain : 782534.203 seconds. 

2. Time used to verify domain P'^': 8797.600 seconds. 

3. Total time: 791331.803 seconds. 

This completes the proof of the problem of spherical distribution of 5 points. 
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A Maple code 

A.l Module Interval Arithmetic 



# IntervalArithmetic : a Maple package used for interval computations. 

# Srevision: 1.0.3.4$ 
IntervalArithmetic := module () 

export ulp, fulb, rfulbO, rfulb, 

'Evalr/add', 'Evalr/multiply' , ' Evalr /power ' , 
'Evalr/sin' , 'Evalr/cos', 'Evalr/tan' , 'Evalr/cot', 
'Evalr/arcsin' , 'Evalr/arccos' , 'Evalr/arctan' , 'Evalr/arccot ' , 
'Evalr/exp', 'Evalr/ln' , ' Evalr /powexp ' , 
' Evalr/shake ' , Evalr, 

isdef , vert, inthull, intwidth, maxwidthdim, intsbdv, sortpos, contain: 
local init: 

option package, load = init: 

############################################################################## 

# initialization 

############################################################################## 

# the initialization procedure 
init := procO 

global 'type/interval', ' type/int_ext ' , 'type/rat_ext ' , 
'type/cons', 'type/consnum' , 'type/ratpar' , 
'type/exp_user' , 'type/In', 

' convert/f t2rat ' , truncatenegativepart , pi, p2 : 

# defining the datatype of interval 
'type/interval' := proc( inv ) 

if type( inv, list( rat_ext ) ) and nops( inv ) = 2 
and inv[ 1 ] <= inv[ 2 ] then 
true 

elif inv = [ ] then 

true 
else 

false 

f i 
end: 

# extended datatype of integer numbers 
'type/int_ext ' := proc( x ) 

evalb( type( x, integer ) or x = -infinity or x = infinity ) 
end: 

# extended datatype of rational numbers 
'type/rat_ext ' := proc( x ) 

evalb( type( x, rational ) or x = -infinity or x = infinity ) 
end: 

# datatype that containing gamma, Catalan or Pi 
'type/cons' := proc( x ) 

member ( x, { gamma, Catalan, Pi } ) 
end: 

# datatype of constant numbers 
'type/consnum' := proc( x ) 

local c, r: 
global constants : 
c := constants: 

constants := gamma, Catalan, Pi, infinity: 
r := type( x, constant ): 
constants := c: 
r 

end: 

# datatype of functions whose parameters are rational numbers 
'type/ratpar' := proc( f, tp::Or( set ( type ), type ) ) 
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type( f , tp ) and type( [ op( f ) ] , list( rational ) ) 
end: 

# datatype of natural exponent 
'type/exp_user ' := proc( f ) 

evalb( op( 0, f ) = 'exp' ) 
end: 

# datatype of natural logarithm 
'type/In' := proc( f ) 

evalb( op( 0, f ) = 'In' ) 
end: 

# converting from float point numbers (ranges) to rational numbers (intervals) 
' convert/f t2rat ' := proc( inv ) 

local t: 

if type( inv, float ) then 

convert ( inv, rational, exact ) 
elif type( inv, rat_ext ) then 

inv 

elif type( inv, list( Or( float, rat_ext ) ) ) and nops( inv ) = 2 then 
t : = map ( procname , inv ) : 
if t[ 1 ] <= t[ 2 ] then t 
else [ ] 
fi 

elif type( inv, list ) or type( inv. Matrix ) then 

map( procname, inv ) 
else 

error "invalid argument: 7,1", inv 

fi 
end: 

# If truncatenegativepart is set to ture, 

# then in interval power arithmetic, 

# negative parts of base intervals are truncated. 

# If truncatenegativepart is set to false, 

# then an error is generated when encountering 

# a interval power arithmetic whose base contains a negative part, 
truncatenegativepart := false: 

# rational lower bound of Pi 
pi := rfulbO( Pi, 'r' , '1' ) : 

# rational upper bound of Pi 
p2 := rfulbO( Pi, 'r' , 'u' ): 

end: 



# the error of a float point number 
ulp := proc( X ) 

if type( X, float ) then 

if X = Float ( infinity ) or x = -Float ( infinity ) then 
else Float( 1, length( op( 1, x ) )+op( 2, x )-Digits ) 
fi 

elif type( x, list ) or type( x. Matrix ) then map( procname, x ) 
else error "invalid argument: '/,!", x 
fi 
end: 

# the float upper or lower bounds of a float number 
fulb := proc( x, ul :: -[identical ( u ), identical ( 1 )}■ ) 

local t: 

if type( X, float ) then 

if x = Float ( infinity ) or x = -Float ( infinity ) then 

X 

else 

# not using "eval" may generate an error, 

# may be a bug in maple ( evaluation rule in Float ). 
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if type( ul, identical ( u ) ) then 

eval( x+Float( 1, length( op( 1, x ) )+op( 2, x )-Digits ) ) 
else 

eval( x-Float( 1, length( op( 1, x ) )+op( 2, x )-Digits ) ) 

fi 

fi 

elif type( x, int_ext ) then 

X 

elif type( x, consnum ) then 

procnameC evalf ( x ), args [ 2..-1 ] ) 
elif type( x, list ) or type( x, Matrix ) then 

map( procname, x, args [ 2..-1 ] ) 
else 

error "invalid argument: 7,1", x 

fi 
end: 

# the rationaKor float) upper(or lower) bound of 

# a number that can be computed to any precision in Maple system 
rfulbO := proc( x, rf :: {identical ( r ), identical ( f )}, 

ul :: -[identical ( u ), identicaK 1 )}■ ) 
local p, q, n, t: 
if type( X, int_ext ) then 
x 

elif type( x, rational ) then 
if rf = 'r' then x 
else fulb( x, ul ) 
f i 

elif type( x, Or( float, cons, 

ratparC { trig, arctrig, exp_user. In } ) ) ) 

or ( op( 0, X ) in { 'exp', 'In', 'arctan', 'arccot' }■ 

and op( x ) = infinity ) 

or ( op( 0, X ) in { 'exp', 'arctan', 'arccot' } 
and op( x ) = -infinity ) then 

if rf = 'f then fulb( x, ul ) 

else convert ( fulb( x, ul ), ft2rat ) 

fi 

elif type( x, ratpar( '"' ) ) then 
q :=op( 1, x): n :=op(2, x): 
if q < then # e.g. ( -16 )"( 1/3 ) 
if type( n, integer ) then 

procnameC q~n, args [ 2..-1 ] ) 
elif type( numer( n ), even ) then 
procnameC C-q)~n, args [ 2..-1 ] ) 
elif typeC numerC n ), odd ) and typeC denomC n ), odd ) then 

procnameC -C -q )~n, args [ 2..-1 ] ) 
else 

error "negative radicand" 

fi 
else 

if rf = 'f then fulbC x, ul ) 

else convert C fulbC x, ul ), ft2rat ) 

fi 

fi 

elif typeC x, '+' ) and 
typeC [ opC X ) ], listC 

OrC rational, cons, ratparC -[trig, arctrig, '"'} ) ) ) ) then 

'+'C opC mapC procname, [ opC x ) ], args [ 2..-1 ] ) ) ) 
elif typeC x, '*' ) and typeC [ opC x ) ] , 
listC OrC rational, ratparC '~' ) ) ) ) then 

if hastypeC removeC type, [ opC x ) ], rational ), negative ) then 

procnameC mapC convert, x, surd ), args [ 2..-1 ] ) 
elif opC select C type, [ opC x ) ], rational ) ) < then 
procnameC mapC procname, x, 'r' , 
'if'C ul = 'u', '1', 'u' ) ), args[ 2..-1 ] ) 
else 

procnameC mapC procname, x, 'r' , 

'if'C ul = 'u', 'u', '1' ) ), args[ 2..-1 ] ) 

fi 
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elif type( x, list ) or type( x, Matrix ) then 

map( procname, x, args [ 2..-1 ] ) 
else 

error "cannot determine the '/,! 7.2 bound of °/,3", 

'if'( rf = 'r', rational, float ), ' if ' ( ul = 'u' , upper, lower ), x 

fi 
end: 

# the rational (or float) upper (or lower) bound of a number 
rfulb := proc( x, rf :: {identical ( r ), identicalC f )}■, 

ul : : -(identicaK u ), identicalC 1 )]■ ) 
local i, t: 

if type( X, consnum ) then 
t := 'Evalr/shake' ( x ): 

if typeC t, 'interval' ) then t := op( ' if ' ( ul = '1', 1, 2 ), t ) fi: 

if rf = 'f then t := rfulbO( t, 'f, ul ) fi: 

t 

elif typeC x, list ) or type( x. Matrix ) then 

map( procname, x ) 
else error "invalid argument: 7,1", x 
fi: 
end: 

############################################################################## 

# Following procedures define interval arithmetic 

############################################################################## 

# interval addition 
'Evalr/add' := proc( L::list ) 

local n, i, lb, ub: 

if member ( [ ] , L) then return [ ] fi: 
n : = nops ( L ) : 
lb := 0: ub := 0: 
for i to n do 

if typeC L[ i ], 'interval' ) then 

lb := lb + L[ i, 1 ] : ub := ub + L[ i, 2 ] : 
elif type( L[ i ], rat_ext ) then 

lb := lb + L[ i ] : ub := ub + L[ i ] : 
else 

error "imrecognized argument: 7,1", L[ i ] 

fi 
od: 

# Large integer may occur after computation with rational numbers, 

# we adjust intervals outside, expressing denominators of interval 

# vertexes with integers whose digits is not greater then Digits, 
lb := convert ( rfulbO( lb, 'f, 1 ), ft2rat ): 

ub := convert ( rfulbO( ub, 'f, u ), ft2rat ): 
return [lb, ub ] 
end: 

# arithmetic multiplication 
'Evalr /multiply' := proc( L::list ) 

local n, i, lb, ub, a, b, tmp: 

if member ( [ ] , L) then return [ ] fi: 

n : = nops ( L ) : 

lb := 1: ub := 1: 

for i to n do 

if type( L[ i ], 'interval' ) then 
a:=L[i, l]:b:=L[i, 2]: 
if a >= then 

if lb >= then lb := lb*a: ub := ub*b: 
elif ub <= then lb := lb*b: ub := ub*a: 
else lb := lb*b: ub := ub*b: 
fi: 

elif b <= then 

if lb >= then tmp := lb: lb := ub*a: ub := tmp*b: 
elif ub <= then tmp := lb: lb := ub*b: ub := tmp*a: 
else tmp := lb: lb := ub*a: ub := tmp*a: 
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fi: 

else # a < < b 

if lb >= then lb := ub*a: ub := ub*b: 
elif ub <= then ub := lb*a: lb := lb*b: 
else # lb < < ub 
tmp : = lb : 

if ( lb*b ) <= ( ub*a ) # lb*b and ub*a will not be 0*infinity 

then lb := lb*b: 
else lb := ub*a: 
fi: 

if ( tmp*a ) <= ( ub*b ) 

then ub : = ub*b : 
else ub := tmp*a: 
fi: 
fi: 
fi: 

elif type( L[ i ], rational ) then 

if ( L [ i ] ) >= then lb : = lb*L [ i ] : ub : = ub*L [ i ] : 

else tmp := lb: lb := ub*L[ i ] : ub := tmp*L[ i ] : 

fi: 

else error "unrecognized argument: L[ i ]: 

f i: 
od: 

lb := convert ( rfulbO( lb, 'f, '1' ), ft2rat ): 
ub := convert ( rfulbO( ub, 'f, 'u' ), ft2rat ): 
return [lb, ub ] 
end: 



# interval power 

' Evalr /power ' := proc( a, n: :rational ) 
local lb, ub, tmp: 
if n = then return fi: 
if a = [ ] then return [ ] fi: 
if type( a, 'interval' ) then 
lb := a[ 1 ] : ub := a[ 2 ] : 
if n > then 
if lb >= then 

lb := lb~n: ub := ub~n: 
elif ub < then 

if type( numer( n ), odd ) then 
if type( denom( n ) , odd ) then 

lb := lb"n: ub := ub~n: 
else # e.g. sqrt ( [ -2, -1 ] ) 

error "negative radicand" 
fi: 

else tmp := lb: lb := ub~n: ub := tmp~n: 
fi: 

else # lb < =< ub 

if type ( numer ( n ) , odd ) then 
if type( denom( n ) , odd ) then 

lb := lb~n: ub := ub'n: 
else # e.g. sqrt( [ -2, 1 ] ) 
if truncatenegativepart then 

lb := 0: ub := ub~n: 
else 

error "negative radicand" 

f i 
fi: 

else # e.g. ( [ -2, 1 ] )~2 

if ( -lb ) <= ub then ub := ub"n: 

else ub := ( -lb )~n: 

fi: 

lb := 0: 
fi: 
fi: 

else # n < 
if lb = then 
if ub = then 

lb := -infinity: ub := infinity: 
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else 

lb := ub~n: ub := infinity: 
fi: 

elif lb > then tmp := lb: lb := ub"n: ub := tmp'n: 
elif ub <= then 

if type( numer( n ), odd ) then 
if type( denom( n ) , odd ) then 

if ub=0 then ub := lb~n: lb := -infinity: 
else tmp := lb: lb := ub~n: ub := tmp~n: 
fi: 

else # e.g. l/sqrt( [ -2,-1 ] ) 

error "negative radicand" 
fi: 
else 

if ub = then lb := lb~n: ub := infinity: 
else lb := lb~n: ub := ub~n: 
fi: 
fi: 

else # lb < < ub 

if type( numer( n ), odd ) then # e.g. l/sqrt( [ -2, 1 ] ) 
if truncatenegativepart then 
lb := ub~n: ub := infinity: 
else 

error "negative radicand" 

fi 

else # e.g. l/( [ -2, 1 ] )'2. 
if ( -lb ) <= ub then lb := ub~n: 
else lb := ( -lb )-n: 
fi: 

ub := infinity: 
fi: 
fi: 
fi: 

elif type( a, rational ) then 

return [ rfulbO( a~n, 'r', '1' ), rfulbOC a"n, 'r' , 'u' ) ] 
else error "invalid argument: 7,1", a 
fi: 

if type( lb, Not( rational ) ) then lb := rfulbO( lb, 'r' , '1' ) 
else lb := convert ( rfulbO( lb, 'f, '1' ), ft2rat ) 
f i: 

if type( ub, Not( rational ) ) then ub := rfulbO( ub, 'r' , 'u' ) 

else ub := convert ( rfulbO( ub, 'f, 'u' ), ft2rat ) 

fi: 

return [lb, ub ] 
end: 



# interval sine 
'Evalr/sin' := proc( a ) 

local lb, ub, tmp: 

global pi, p2: 

if a = [ ] then return [ ] fi: 
if type( a, 'interval' ) then 
lb := a[ 1 ] : ub := a[ 2 ] : 
if ub <= then 

tmp := procnameC [ -ub, -lb ] ): 
return [ -tmp [ 2 ] , -tmp [ 1 ] ] 
f i: 

lb := lb/'if'( lb >0, p2, pi ) : ub := ub/'if'( ub >0, pi, p2 ): 
tmp := floor ( lb/2 ): lb := lb-2*tmp: ub := ub-2*tmp: 
if ub >= ( lb+2 ) then lb := -1: ub := 1: 
else 

if ( lb-1/2 ) < then 
if ( ub-1/2 ) < then 

# l/2*Pi > lb*Pi > lb*pl > 

# ( ub-1/2 )*Pi < ( ub-1/2 )*pl < ==> 

# < ub*Pi < l/2*Pi+( ub-1/2 )*pl < Pi/2 
lb := rfulbOC sin( lb*pl ), 'r' , '1' ): 

ub := rfulbOC cos( ( 1/2-ub )*pl ), 'r' , 'u' ): 
elif ( ub-l+lb ) < then 
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'1' ): ub 



# l/2*Pi > lb*Pi > lb*pl > 
lb := rfulbO( sinC lb*pl ), ' 

elif ( ub-1 ) < then 

# ( ub-1 )*Pi < ( ub-1 )*pl < ==> 

# Pi/2 < ub*Pi < Pi+( ub-1 )*pl < Pi 
lb := rfulbOC sin( ( 1-ub )*pl ), 'r 

elif ( ub-3/2 ) < then 

# ( ub-3/2 )*Pi < ( ub-3/2 )*pl < ==> 

# Pi <= ub*Pi < 3/2*Pi+( ub-3/2 )*pl < 3/2*Pi 
lb := -rfulbOC cos( ( 3/2-ub )*pl ), 'r', 'u' 

else 

lb := -1: ub := 1: 
fi: 

elif ( lb-1 ) < then 
if ( ub-3/2 ) < then 

# ( ub-3/2 )*Pi < ( ub-3/2 )*pl <0 ==> 

# Pi/2 <= ub*Pi < 3/2*Pi+( ub-3/2 )*pl 

# ( lb-1/2 )*Pi >= ( lb-1/2 )*pl >=0 = 



1: 



) : ub : = 1 : 



) : ub 



3/2*Pi 



# Pi > lb*Pi >= Pi/2+( lb-1/2 )*pl >= Pi/2 
tmp : = lb : 

lb := -rfulbOC cos( ( 3/2-ub )*pl ) 



'r ' , 



ub := rfulbOC cos( ( tmp-1/2 )*pl ). 
elif ( ub-3+lb ) < then 

# ( lb-1/2 )*Pi >= ( lb-1/2 )*pl >=0 ==> 

# Pi > lb*Pi >= Pi/2+( lb-1/2 )*pl >= Pi/2 

ub := rfulbOC cos( ( lb-1/2 )*pl ), 'r' , 'u' ) : lb := -1: 
elif ( ub-5/2 ) < then 

# ( ub-5/2 )*Pi < ( ub-5/2 )*pl < ==> 

# 2*Pi <= ub*Pi < 5/2*Pi+( ub-5/2 )*pl < 5/2*Pi 

lb := -1: ub := rfulbO( cos( ( 5/2-ub )*pl ), 'r', 'u' ): 
else lb := -1: ub := 1: 
fi: 

elif ( lb-3/2 ) < then 
if ( ub-3/2 ) < then 

# ( ub-3/2 )*Pi < ( ub-3/2 )*pl < ==> 

# Pi < ub*Pi < 3/2*Pi+( ub-3/2 )*pl < 3/2*Pi 

# ( lb-1 )*Pi >= ( lb-1 )*pl >= ==> 

# 3/2*Pi > lb*Pi >= Pi+( lb-1 )*pl >= Pi 
tmp : = lb : 

lb := -rfulbOC cos( ( 3/2-ub )*pl ), 'r', 
ub := -rfulbOC sin( ( tmp-1 )*pl ), 'r' , 
elif ( ub-3+lb ) < then 

# ( lb-1 )*Pi >= ( lb-1 )*pl >=0 ==> 

# 3/2*Pi > lb*Pi >= Pi+( lb-1 )*pl >= Pi 
ub := -rfulbOC sin( ( lb-1 )*pl ), 'r', ' 

elif ( ub-5/2 ) < then 

# ( ub-5/2 )*Pi < ( ub-5/2 )*pl < ==> 

# 3/2*Pi <= ub*Pi < 5/2*Pi+( ub-5/2 )*pl < 5/2*Pi 
lb := -1: ub := rfulbO( cos( ( 5/2-ub )*pl ), 'r', 

else lb := -1: ub := 1: 
fi: 

elif ( ub-5/2 ) < then 

# ( lb-3/2 )*Pi >= ( lb-3/2 )*pl >=0 ==> 

# 2*Pi > lb*Pi >= 3/2*Pi+( lb-3/2 )*pl >= 3/2*Pi 

=> 



' ) 
): 



): lb := 



# ( ub-5/2 )*Pi < ( ub-5/2 )*pl < 

# 3/2*Pi <= ub*Pi < 5/2*Pi+( ub-5/2 )*pl 
lb := -rfulbOC cos( ( lb-3/2 )*pl ), 'r' 
ub := rfulbOC cos( ( 5/2-ub )*pl ), 'r' , 

elif ( ub-5+lb ) < then 

# ( lb-3/2 )*Pi >= ( lb-3/2 )*pl >=0 ==> 

# 2*Pi > lb*Pi >= 3/2*Pi+( lb-3/2 )*pl >= 

), 'r' 



lb := -rfulbOC cos( ( lb-3/2 )*pl 
elif ( ub-7/2 ) < then 

# ( ub-7/2 )*Pi < ( ub-7/2 )*pl < 

# 3*Pi <= ub*Pi < 7/2*Pi+( ub-7/2 
lb := -rfulbOC cos( ( 7/2-ub )*pl 

else lb := -1: ub := 1: 
fi: 
fi: 



< 5/2*Pi 

'u' ): 
'u' ): 



3/2*Pi 
'u' ) : ub 



==> 
)*pl < 
), 'rV 



7/2*Pi 
'u' ): 



ub 



1: 



1: 
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elif type( a, rational ) then 

return [ rfulbO( 'sin'( a ), 'r' , '1' ), rfulbO( 'sin'( a ), 'r', 'u' ) ] 
else error "invalid argument: 7,1", a 
fi: 

return [lb, ub ] 
end: 

# interval cosine 
'Evalr/cos' := proc( a ) 

global pi, p2: 

if a = [ ] then return [ ] fi: 

if type( a, 'interval' ) then 

'Evalr/sin'( [ -a[ 2 ] + pl/2, -a[ 1 ] + p2/2 ] ) 

elif type( a, rational ) then 

# using 'cos' to suppress automatic simplification 

[ rfulbO( 'cos'( a ), 'r', '1' ), rfulbO( 'cos'( a ), 'r', 'u' ) ] 

else 

error "invalid argument: '/,!", a 
fi: 
end: 

# interval tangent 
'Evalr/tan' := proc( a ) 

local lb, ub, tmp: 
global pi , p2 : 

if a = [ ] then return [ ] fi: 
if type( a, 'interval' ) then 
lb := a[ 1 ] : ub := a[ 2 ] : 

lb := lb/'if'( lb > 0, p2, pi ) : ub := ub/'if'( ub > 0, pi, p2 ): 

tmp := floor ( lb ) : lb := lb - tmp: ub := ub - tmp: 

if ub >= ( lb+1 ) then lb := -infinity: ub := infinity: 

else 

if ( lb-1/2 ) < then 
if ( ub-1/2 ) < then 

# 0<lb*pl<lb*Pi<a [1] -tmp*Pi<a [2] -tmp*Pi<ub*Pi<ub*p2 

lb := rfulbOC 'tan'( a[ 1 ] ), 'r' , '1' ): 

ub := rfulbO( 'tan'( a[ 2 ] ), 'r' , 'u' ): 
else 

lb := -infinity: ub := infinity: 
fi: 

elif ( ub-3/2 ) < then 

# 0<lb*pl<lb*Pi<a [1] -tmp*Pi<a [2] -tmp*Pi<ub*Pi<ub*p2 

lb := rfulbOC 'tan'( a[ 1 ] ), 'r', '1' ): 

ub := rfulbOC 'tan'( a[ 2 ] ), 'r', 'u' ): 
else 

lb := -infinity: ub := infinity: 
fi: 
fi: 

elif type( a, rational ) then 

return [ rfulbOC 'tan'C a ), 'r' , '1' ), rfulbOC 'tan'C a ), 'r', 'u' ) ] 
else 

error "invalid argument: 7,1", a 
fi: 

return [lb, ub ] 
end: 

# interval cotangent 
'Evalr/cot' := procC a ) 

global pi , p2 : 

if a = [ ] then return [ ] fi: 
if typeC a, 'interval' ) then 

'Evalr/tan' C [ -a[ 2 ] + pl/2, -a[ 1 ] + p2/2 ] ) 
elif typeC a, rational ) then 

[ rfulbOC 'cot'C a ), 'r', '1' ), rfulbOC 'cot'C a ), 'r', 'u' ) ] 
else 

error "invalid argument: 7ol", a 
fi: 
end: 
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# interval arc sine 
'Evalr/arcsin' := proc( inv ) 

local lb, ub: 

global pi , p2 : 

if inv = [ ] then return [ ] fi: 
if type( inv, rational ) then 
if inv >= -1 and inv <= 1 then 

[ rfulbOC 'arcsin'C inv ), 'r', '1' ), 
rfulbOC 'arcsin'C inv ), 'r' , 'u' ) ] 
else 

error "invalid argument: 7,1", inv: 

fi 

elif type( inv, 'interval' ) then 
lb : = inv [ 1 ] : ub : = inv [ 2 ] : 

if not ( type( lb, rational ) and type( ub, rational ) ) then 

error "invalid argument: 7,1", inv 
fi: 

if lb > 1 or ub < -1 then 

error "invalid argument: 7«1", inv: 
elif lb < -1 and ub > 1 then 

lb := -p2: 

ub := p2: 
elif lb < -1 then 

lb := -p2 : ub := rfulbOC 'arcsin'C ub ), 'r', 'u' ): 
elif ub > 1 then 

lb := rfulbOC 'arcsin'C lb ), 'r' , '1' ) : ub := p2: 
else 

lb := rfulbOC 'arcsin'C lb ), 'r' , '1' ): 
ub := rfulbOC 'arcsin'C ub ), 'r' , 'u' ): 
f i: 

[ lb, ub ] 
else 

error "invalid argument: 7.1", inv 

fi 
end: 



# interval arc cosine 
'Evalr/arccos' := procC inv ) 

local tmp: 

global pi , p2 : 

if inv = [ ] then return [ ] fi: 
if typeC inv, rational ) then 
if inv >= -1 and inv <= 1 then 

[ rfulbOC 'arccos'C inv ), 'r' , '1' ), 
rfulbOC 'arccos'C inv ), 'r' , 'u' ) ] 
else 

error "invalid argument: 7.1", inv: 
fi: 

elif typeC inv, 'interval' ) then 
tmp := 'Evalr/arcsin' C inv ): 
[ pl/2-tmp[ 2 ], p2/2-tmp[ 1 ] ] 

else 

error "invalid argument: 7.1", inv: 
fi: 
end: 



# interval arc tangent 
'Evalr/arctan' := procC inv ) 

if inv = [ ] then return [ ] fi: 

if typeC inv,rat_ext ) then 

[ rfulbOC 'arctan'C inv ), 'r' , '1' ), 
rfulbOC 'arctan'C inv ), 'r', 'u' ) ] 

elif typeC inv, 'interval' ) then 

[ rfulbOC 'arctan'C inv[ 1 ] ), 'r' , '1' ), 
rfulbOC 'arctan'C inv [ 2 ] ), 'r', 'u' ) ] 

else 

error "invalid argument: 7.1", inv: 

fi 
end: 



25 



# interval arc cotangent 
'Evalr/arccot ' := proc( inv ) 

local tmp: 

global pi, p2: 

if inv = [ ] then return [ ] fi: 

if type( inv, rat_ext ) then 

[ rfulbOC 'arccot'C inv ), 'r' , '1' ), 
rfulbOC 'arccot'C inv ), 'r', 'u' ) ] 

elif type( inv, 'interval' ) then 
tmp := 'Evalr/arctan' ( inv ): 
[ pl/2-tmp[ 2 ], p2/2-tmp[ 1 ] ] 

else 

error "invalid argument: '/,!", inv: 

fi 
end: 



# interval exponent 
'Evalr/exp' := proc( ) 
local a, inv: 

if has( [args] , { [ ] } ) then return [ ] fi: 
if nargs = 1 then 
inv : = args [ 1 ] : 
if type( inv, rat_ext ) then 

[ rfulbOC exp( inv ), 'r' , '1' ), rfulbOC exp( inv ), 'r', 'u' ) ] 
elif inv = [ ] then 
[ ] 

elif type( inv, 'interval' ) then 

[ rfulbOC exp( inv[ 1 ] ), 'r' , '1' ), 

rfulbOC expC inv[ 2 ] ) , 'r ' , 'u' ) ] 
else 

error "invalid argument: 7,1", inv: 

f i 

elif typeC args [ 1 ] , rational ) and args [ 1 ] > and args [ 1 ] <> 1 then 
a : = args [ 1 ] : 
inv : = args [ 2 ] : 
if typeC inv, rat_ext ) then 

[ rfulbOC a" inv, 'r' , '1' ), rfulbOC a" inv, 'r', 'u' ) ] 
elif inv = [ ] then 

[ ] 

elif a > 1 and typeC inv, 'interval' ) then 

[ rfulbOC a'C inv[ 1 ] ), 'r', '1' ), 

rfulbOC a'C inv[ 2 ] ) , 'r ' , 'u' ) ] 
elif a < 1 and typeC inv, 'interval' ) then 

[ rfulbOC a'C inv[ 2 ] ), 'r', '1' ), 

rfulbOC a'C inv[ 1 ] ) , 'r ' , 'u' ) ] 
else 

error "invalid argument: 7,1, 'Li", a, inv: 

fi 
else 

error "invalid argument: 7,0", args: 

fi 
end: 



# interval logarithm 
'Evalr/ln' := procC ) 
local a, inv: 

if hasC [args] ,{[]}) then return [ ] fi: 
if nargs = 1 then 
inv : = args [ 1 ] : 

if typeC inv, rat_ext ) and inv > then 

# e.g. lnC4) is automatically simplified to 2*lnC2) . 

[ rfulbOC 'In'C inv ), 'r', '1' ), rfulbOC 'In'C inv ), 'r' , 'u' ) ] 
elif inv = [ ] then 

[ ] 

elif typeC inv, 'interval' ) and inv[ 1 ] >= then 

[ 'if'C inv[ 1 ] = 0, -infinity, rfulbOC 'In'C inv[ 1 ] ), 'r', '1' ) ), 

rfulbOC 'In'C inv[ 2 ] ) , 'r ' , 'u' ) ] 
else 
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error "invalid argument: 7,1", inv: 

fi 

elif type( args [ 1 ] , rational ) and args [ 1 ] > and args [ 1 ] <> 1 then 
a : = args [ 1 ] : 
inv : = args [ 2 ] : 

if type( inv, rat_ext ) and inv > then 
[ rfulbOC 'ln'( inv ), 'r', '1' )/ 

rfulbO( 'ln'( a ), 'r' , ' if ' ( inv > 1, 'u' , '1' ) ), 
rfulbOC 'ln'( inv ), 'r', 'u' )/ 

rfulbOC 'ln'( a ), 'r' , ' if ' ( inv > 1, '1', 'u' ) ) ] 
elif inv = [ ] then 
[ ] 

elif a > 1 and type( inv, 'interval' ) and inv[ 1 ] >= then 
[ 'if'( inv[ 1 ] = 0, -infinity, 
rfulbOC 'ln'( inv[ 1 ] ), 'r', '1' )/ 

rfulbOC 'ln'( a ), 'r' , ' if ' ( inv[ 1 ] > 1, 'u' , '1' ) ) ), 
rfulbOC 'In'C inv[ 2 ] ), 'r', 'u' )/ 

rfulbOC 'In'C a ), 'r' , ' if ' C inv[ 2 ] > 1, '1', 'u' ) ) ] 
elif a < 1 and typeC inv, 'interval' ) and inv[ 1 ] >= then 
[ rfulbOC 'In'C inv[ 2 ] ), 'r', '1' )/ 

rfulbOC 'In'C a ), 'r' , ' if ' C inv[ 2 ] > 1, 'u' , '1' ) ), 
'if'C inv[ 1 ] = 0, infinity, 
rfulbOC 'In'C inv[ 1 ] ), 'r', 'u' )/ 

rfulbOC 'In'C a ), 'r' , 'if'C inv[ 1 ] > 1, '1', 'u' ) ) ) ] 
else 

error "invalid argument: 7,1, 7.2", a, inv: 

fi 
else 

error "invalid argument: 7,2", args: 

fi 
end: 



# interval power exponent 

'Evalr/powexp' := procC a::interval, x::interval ) 
ifa= [] orx= [] then 

[ ] 

elif a = and x[ 1 ]*x[ 2 ] > then 

[ 0, ] 
elif a = 1 then 

[ 1, 1 ] 

# elif a[ 1 ] >= and 'if'C a[ 1 ] = 0, x[ 1 ]*x[ 2 ] >= 0, true ) then 
elif a[ 1 ] >= then 

'Evalr/exp'C 'Evalr/multiply ' C [ x, 'Evalr/ln' C a ) ] ) ) 
else 

error "invalid argument in 7il: 7.2, 7o3", procname, a, x: 

fi 
end: 



# to get an interval containing the number 
' Evalr/shake ' := procC expr ) 
local tm, tp: 

if Digits < 20 then Digits := 2*Digits: fi: 

tm := [ opC expr ) ] : 

tp := opC 0, expr ) : 

if typeC expr, rat_ext ) then 

[ expr, expr ] 
elif typeC expr, QrC float, cons ) ) then 

[ rfulbOC expr, 'r' , '1' ), rfulbOC expr, 'r' , 'u' ) ] 
elif tp = '+' then 

'Evalr/add'C mapC procname, tm ) ) 
elif tp = '*' then 

'Evalr/multiply ' C mapC procname, tm ) ) 
elif tp = ' ~ ' then 

if typeC tm[ 2 ], rational ) then 

'Evalr/power ' C procnameC tm[ 1 ] ), tm[ 2 ] ) 

elif typeC tm[ 1 ], rational ) then 

'Evalr/exp'C tm[ 1 ], procnameC tm[ 2 ] ) ) 

else 

'Evalr/powexp' C procnameC tm[ 1 ] ), procnameC tm[ 2 ] ) ) 
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fi 

elif tp in ■[' sin' , ' cos ' , ' tan' , ' cot ' , 

' arcsin' , ' arccos ' , ' arctan' , ' arccot ' , ' exp' , ' In' }■ then 

thismodule [ cat ( 'Evalr/',tp ) ]( procname( tm[ 1 ] ) ) 
elif type( expr, Qr( list, Matrix ) ) then 

map( procname, expr ) 
else 

error "invalid argument: /ll", expr 

fi 
end: 



# interval arithmetic 

Evalr := proc( expr, invls::Or( list, set ) ) 
local i, vars , t, tm, tp, tmp: 
option system, remember: 

vars := select ( type, indets( expr ), name ): 
if vars = { } then 

if type( expr, consnum ) then 
return ' Evalr/ shake ' ( expr ) 
elif has ( expr ,-[[]}) then 

return [ ] 
elif type( expr, 'interval' ) then 
return expr 

elif type( expr, list ) and nops( expr ) = 2 and 

type ( expr [ 1 ] , consnum ) and type ( expr [ 2 ] , consnum ) then 

return [ ' Evalr/ shake ' ( expr [ 1 ] ) [ 1 ] , 

'Evalr/shake' ( expr[ 2 ] )[ 2 ] ] 
elif type( expr, Or( list. Matrix ) ) then 

return map( procname, expr, args [ 2.. -1 ] ) 
elif not member ( op( 0, expr ), 
{'+', '*', '-', 'sin', 'cos', 'tan', 'cot', 
'arcsin', 'arccos', 'arctan', 'arccot', 'exp', 'In'} ) then 

error "invalid argument: 7,1", expr: 
fi: 
else 

if nargs = 1 then 

return procname ( expr, map( t -> t = [ -infinity, infinity ] , vars ) ) 
elif member ( [ ], subs( invls, vars ) ) then 
return [ ] 

elif vars minus ■[ op( map( Ihs , invls ) )}■<>{} then 
return procname ( expr, 

{ op( invls ) } union map( t -> t = [ -infinity, infinity ] , 
vars minus { op( map( Ihs, invls ) ) } ) ) 
elif not type( subs( invls, vars ), set ( 'interval' ) ) then 
tmp : = [ ] : 

t := [ op( map( t -> t = subs( invls, t ), vars ) ) ] : 
for i to nops( t ) do 

if type( rhs( t[ i ] ), 'interval' ) then 

tmp : = [ op ( tmp ) , t [ i ] ] : 
elif type( rhs( t[ i ] ), list ) and nops( rhs( t[i] ) ) = 2 and 
type( rhs( t[i] )[1], consnum ) and 
type ( rhs (t[i] )[2], consnum ) then 

tmp : = [ op ( tmp ) , Ihs ( t [ i ] ) = 

[ 'Evalr/shake' ( rhs( t[i] )[1] )[1], 

'Evalr/shake' ( rhs( t[ i ] )[ 2 ] )[ 2 ] ] ]: 
elif type( rhs( t[ i ] ), consnum ) then 

tmp : = [ op ( tmp ) , ' Evalr/ shake ' ( rhs (t[i] ) ) ]: 
else 

error "invalid evaluation range: 7,1", rhs( t[ i ] ): 
fi: 
od: 

return procname ( expr, tmp ) 
elif type( expr, name ) then 
return subs( invls, expr ) 
fi: 
fi: 



tm := [ op( expr ) ] : 
tp := op( 0, expr ) : 
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if tp = '+' then 

'Evalr/add'( map( procname, tm, args [ 2..-1 ] ) ) 
elif tp = '*• then 

'Evalr/multiply ' ( map( procname, tm, args [ 2..-1 ] ) ) 
elif tp = ' ~ ' then 

if type( tm[ 2 ], rational ) then 

'Evalr/power ' ( procname( tm[ 1 ], args [ 2..-1 ] ), tm[ 2 ] ) 

elif type( tm[ 1 ] , rational ) then 

'Evalr/exp'C tm[ 1 ], procname( tm[ 2 ], args [ 2..-1 ] ) ) 

else 

'Evalr/powexp' ( procname( tm[ 1 ], args [ 2..-1 ] ), 
procname ( tm[ 2 ], args [ 2..-1 ] ) ) 

fi 

elif tp in {'sin', 'cos', 'tan', 'cot', 

'arcsin', 'arccos', 'arctan', 'arccot', 'exp', 'In'}- then 

thismodule[ cat ( 'Evalr/',tp ) ]( procname( tm[ 1 ], args [ 2..-1 ] ) ) 

elif type( expr, Or( list. Matrix ) ) then 
map( procname, expr, args [ 2.. -1 ] ) 

else 

error "invalid argument: 7,1 ", expr 

fi 
end: 

############################################################################## 

# Auxiliary procedures 

############################################################################## 

# list sorting 

sortpos := proc( Is:: list ) 
local V, i, j, n, tmp: 
n := nops( Is ) : 

V := [ seq( [ Is [ i], i], i=l..n) ]: 
for i from 2 to n do 

if ( v[ i ] [ 1 ] ) < ( v[ i-1 ] [ 1 ] ) then 

tmp : = V [ i ] : 

v[ i ] := v[ i-1 ] : 

for j from i-2 by -1 to 1 while ( tmpC l])<(v[j][l])do 

v[ j+1 ] := v[ j ] : 
od: 

if j = 1 and ( tmp [ 1] )<(v[j][l] ) then v[ j ] := tmp: 
else v[ j+1 ] := tmp: 
fi: 
fi: 
od: 

return map2( op, 2, v ) : 
end: 



# The "isdef" procedure determines whether or not a interval matrix is 

# positive definite 

# input: interval matrix and options, can receive 1 to 3 parameters 

# interval matrix can have following forms: 
## 1 . ( [ a, b ] )_{nn}, one parameter 

## 2. [ A, B ] , two parameters 

# option query=posdef /negdef : 

## determine positive (default) /negative def initeness . 

# option query=possemidef /negsemidef : 

## determine positive /negative semidef initeness . 

# option query=nonposdef /nonnegdef : 

## determine nonpositive/nonnegative def initeness . 

# option query=nonpossemidef /nonnegsemidef : 

## determine nonpositive/nonnegative semidef initeness . 

# option query=nondef : 

## determine nondef initeness . 

# option method=vertex/eigenvalue: 

## use different methods to determine positive /negative def initeness : 

## the method using vertex matrices(a sufficient and necessary condition, 

## the default method), and the method using eigenvalues (a sufficient condition) 

## in which case the "false" returned just indicates an unsuccessful test. 

# output: true/false 
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# Note: in order to reduce the time of datatype checking, 

# the entries in matrices should have the same datetype. 
isdef := procO 

local A, B, n, _query, _method, i, j, k, q, S, v, Ac, Ad, Ace, Ade, temp, x: 
uses Linear Algebra: 

A :=0: B :=0: n :=0: _query := 0: .method := 0: 
for i to nargs do 

if type( args [ i ], 'Matrix' ( square ) ) then 

# more accurate form should be: 

# if type( args [ i ], 'Matrix' ( interval ) ) then 
if type( args [ i ][ 1, 1 ], 'interval' ) then 

A : = map2 ( op , 1 , args [ i ] ) : 
B : = map2 ( op , 2 , args [ i ] ) : 
n := RowDimensionC A ) : 

# more accurate form should be: 

# elif type( args [ i ], 'Matrix' ( constant ) ) then 

# and comparison should be contained. 

elif type( args [ i ][ 1, 1 ], rat_ext ) then 
if A = then 

A : = args [ i ] : 

n := RowDimensionC A ): 
else 

if RowDimensionC args [ i ] ) <> n then 

error "different matrix dimension: 7,1", args [ i ] 
fi: 

B : = args [ i ] : 
fi: 

else error "error matrix: 7,1", args [ i ]: 
fi: 

elif typeC args [ i ], 
identical C query ) = 

■[ identical C posdef ), identical C negdef ), 
identicalC possemidef ), identicalC negsemidef ), 
identical C nonposdef ) , identicalC nonnegdef ) , 
identicalC nonpossemidef ), identicalC nonnegsemidef ), 
identicalC nondef ) } ) then 

_query : = rhs C args [ i ] ) : 
elif typeC args [ i ], identicalC method ) = 
■[ identicalC vertex ), identicalC eigenvalue ) } ) then 

_method : = rhs C args [ i ] ) : 
else error "invalid argument: 7,1", args [ i ]: 
f i: 
od: 

ifA = 0orB = 0orn = then error "invalid arguments: 7oO", args fi: 
if hasC A, infinity ) or hasC B, infinity ) then return false fi: 
if _query = then _query := posdef fi: 
if member C _query, 

{ 'posdef, 'negdef, 'possemidef, 'negsemidef } ) and 
_method = then 

_method := vertex 
fi: 

if member C _query, 

{ nonposdef , nonnegdef , nonpossemidef , nonnegsemidef , nondef } ) and 
_method = 'vertex' then 

error "the method specified cannot be applied" 
fi: 

if _method = 'vertex' then 
q := 'positive_def inite ' ; 

if _query = 'negdef or _query = 'negsemidef then 

temp := A: A := -B: B := -temp: 
fi: 

if _query = 'possemidef or _query = 'negsemidef then 
q := 'positive_semidef inite ' ; 

fi; 

for k from to 2"C n-1 )-l do 
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V := VectorC n, convert( k, base 
S := Matrix ( n, n, ( i, j ) -> 
( ( A[ i, j ]+B[ i, j ] )+ 
( -1 )-( v[ i ]+v[ j ] )*( A[ i, j 
shape = symmetric ) : 
if LA_Main:-IsDef inite( S, 'query' 

return false 
fi: 
od: 

return true 
else 

Ac := Matrix ( n, n, (i, j)->(A[i, 
shape = symmetric ) : 

Ad := Matrix ( n,n, (i, j)->(-A[i, 
shape = symmetric ) 
Ace 



2 ) ): 

i ]-B[ i, j ] ) )/2, 
q ) = false then 



j ]+B[ i, 
j ]+B[ i. 



] )/2, 
j ] )/2, 



realrootC CharacteristicPolynomiaK Ac, x 
Ade := realrootC CharacteristicPolynomiaK Ad, x 
Ade := max( op( map( abs, map( op, Ade ) ) ) ): 
if _query = 'posdef and 

min( op( map2( op, 1, Ace ) ) )-Ade > then 

return true 
elif _query = 'negdef and 

Ace ) ) )+Ade < then 



10 
10 



■( -Digits 
■( -Digits 



' and 
) ) )■ 

' and 



Ade >= then 



'nonposdef ' 
op, 1, Ace 

' nonnegdef ' 
op, 2, Ace 



and 
) ) 

and 
) ) 



)+Ade <= then 



)-Ade >= then 



max( op( map2( op, 2 

return true 
elif _query = 'possemidef 
min( op( map2( op, 1, Ace 

return true 
elif _query = 'negsemidef 
max( op( map2( op, 2, Ace ) ) )+Ade <= then 

return true 
elif _query = 
min( op( map2( 

return true 
elif _query = 
max( op( map2( 

return true 
elif _query = 
min( op( map2( 

return true 
elif _query = 

max( op( map2( op, 2, Ace ) ) 

return true 
elif min( op( map2( op, 1 
max( op( map2( op, 2, Ace 

return true 
else 

return false 
fi: 
fi: 
end: 



'nonpossemidef ' 
op, 1, Ace ) ) 

'nonnegsemidef ' 



and 
)+Ade 

and 
)-Ade 



< then 



> then 



, Ace ) ) )+Ade < 
) ) )-Ade > then 



and 



# get vertexes of a domain 
vert := proc( va::list ) 

local n,i,j,k,vt: 

n : = nops ( va ) : 

vt := []: 

if type( va[ 1 ] , '=' ) then 

for i from to 2~n-l do 

j := VectorC n, convertC i, base, 2 ) ): 

vt := [ op( vt ), [ seq( lhs( va[ k ] ) 

=op( j[ k rhs( va[ k ] ) ), k=l..n ) ] ]: 

od: 

elif type( va[ 1 ], 'interval' ) then 
for i from to 2~n-l do 

j := VectorC n, convertC i, base, 2 ) ): 

vt := [ opC vt ) , [ seqC opC j[ k va[ k ] ), k=l..n ) ] ] 

od: 
fi: 

return vt : 
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end: 



# get the convex hull of a set of intervals 
inthull := proc( invls :: list (interval) ) 

if nops ( invls ) = 1 then 
invls [ 1 ] 

else 

[ op( sortposC map2( op, 1, invls ) )[ 1 ], invls )[ 1 ], 
op( sortposC map2( op, 2, invls ) ) [ -1 ] , invls )[ 2 ] ] 

fi 
end: 



# get the width of an interval 
intwidth := proc( inv ) 

if type( inv, 'interval' ) then 

inv[ 2 ] - inv[ 1 ] 
elif type( inv, rational ) then 


elif type( inv, '=' ) and type( rhs( inv ), 'interval' ) then 

rhs ( inv ) [ 2 ] - rhs ( inv ) [ 1 ] 
else 

error "invalid argument: '/,!", inv 

fi 
end:# 



# get the position of the interval with the maximal width in 

# a list of intervals 
maxwidthdim := proc( invls:: list ) 

op( -1, sortposC map( intwidth, invls ) ) ) 
end: 



# subdivides a domain over an interval or all intervals 
intsbdv := proc( invls:: list, sbdv :: integer ) 

local i, j, dim, s, t, n, cnt , subint : 

n := nops( invls ): 

if sbdv = 1 or has( invls, infinity ) then return [ invls ] fi: 
dim := 0: 

if nargs > 2 and type( args [ 3 ] , integer ) and args [ 3 ] <= n then 

dim : = args [ 3 ] 
elif nargs > 2 then 

error "invalid argument: 7,1", args [ 3 ] 
fi: 

if dim <> then 

if type( invls [ 1 ], '=' ) then 

return [ seq( [ op( invls [ l..dim-l ] ), lhs( invls [ dim ] ) = 
[ rhs( invls [ dim ] )[ 1 ] + 
i/sbdv*intwidth( rhs( invls [ dim ] ) ), 
rhs( invls [ dim ] ) [ 1 ] + 

( i+1 ) /sbdv*intwidth( rhs( invls [ dim ] ) ) ], 
op( invls [ dim+l..n ] ) ], i = 0..sbdv-l ) ]: 
elif type( invls [ 1 ], 'interval' ) then 
return [ seq( [ op( invls [ l..dim-l ] ), 

[ invls [ dim, 1 ] + i/sbdv* intwidth ( invls [ dim ] ), 
invls [ dim, 1 ] + ( i+1 ) /sbdv*intwidth( invls [ dim ] ) ], 
op( invls [ dim+l..n ] ) ], i = 0..sbdv-l ) ]: 

fi: 
fi: 

cnt := sbdv"n-l: 
subint : = [] : 

if type( invls [ 1 ], '=' ) then 
for i from to cnt do 

j := Vector[ row ]( n,convert( i, base, sbdv ) ): 
subint : = [ op ( subint ) , 

[ seq( lhs( invls [ s ] ) = 
[ rhs( invls [ s ] ) [ 1 ] + 

j[ s ]/sbdv*( rhs( invls [ s ] )[ 2 ]-rhs( invls [ s ] )[ 1 ] ), 
rhs ( invls [ s ] ) [ 1 ] + 

( j[ s ]+l )/sbdv*( rhs( invls [ s ] )[ 2 ]- 
rhs ( invls [ s ] ) [ 1 ] ) ] , 
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s=l..n ) ] ]: 

od 

elif type( invls [ 1 ], 'interval' ) then 
for i from to cnt do 

j := Vector[ row ] ( n, convert( i, base, sbdv ) ): 
subint : = [ op ( subint ) , 
[ seq( 

[ invls [ s ] [ 1 ] + 

j[ s ]/sbdv*( invls [ s ][ 2 ] -invls [ s ][ 1 ] ), 
invls [ s ] [ 1 ] + 

( j[ s ]+l )/sbdv*( invls[ s ][ 2 ]-invls[ s ][ 1 ] ) ], 
s=l. .n ) ] ] : 

od 
fi: 

return subint 
end: 

# test whether a rational number or an interval is contained in an interval 
contain := proc( inv::Or( name = interval, interval ), 
p::Or( name = Or( rational, interval ), rational, interval ) ) 
if type( inv, name = 'interval' ) then 

procname ( rhs ( inv ) , p ) 
elif type( p, name = Or( rational, 'interval' ) ) then 

procname ( inv, rhs( p ) ) 
elif type( p, rational ) then 

evalb ( inv [ 1 ] <= p and p <= inv [ 2 ] ) 
elif type( p, 'interval' ) then 

evalb ( inv[ 1 ] <= p[ 1 ] and p[ 1 ] <= inv[ 2 ] and 
inv[ 1 ] <= p[ 2 ] and p[ 2 ] <= inv[ 2 ] ) 
else 

error "invalid arguments: 7.0", args 

fi 
end: 

end module : 
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A. 2 Module fivepoints 



# fivepoints: a Maple package used for solving the problem of spherical 

# distribution of 5 points. 

# Srevision: 1.0$ 
fivepoints := module () 

export spf , cmb, changecoord, eucdis, grad, normal vector , 
ischecked, spchecked, 

howchecked, inwhichdomain, subdivisiondomain, subdivisionprocess : 
local init : 

option package, load = init: 
init := proc( ) 

local p, q, i, j, k, m, nv, temp, x, y, A, B, C, D, E, py, pyd. Is: 
global truncatenegativepart , 

phi, theta, indet , Vt , Vtd, f, F, Fls, Df , H, 

Thetabp, fmax, f maxrat , pyang, Thetap, pyfmax, xi, chi, eta, zeta, bpyva, pyva, 
varsrangel, varsrange2, ckrangelsl, ckrangels2, 
gr, grcv, signdis, Methods: 
uses IntervalArithmetic : 



# Negative parts of base intervals in interval power arithmetic are truncated. 

# e.g. Evalr( op( 6, f ), ckrangelsl [ 1 ] ); 
truncatenegativepart := true: 

if FileTools[ Exists ]( "fivepoints . data" ) = false then 



# spherical coordinates of five points 

Vt := [ [ 1, 0, ] , [ 1, phi[ 1 ] , Pi ] , [ 1, phi[ 2 ], theta[ 2 ] ], 
[ 1, phi[ 3 ], theta[ 3 ] ], [ 1 , phi [ 4 ] , theta[ 4 ] ] ]: 

# Cartesian coordinates corresponding the spherical coordinates 
Vtd := map( changecoord, Vt , s2d ): 

# indeterminates 

indet := [ phi [ 1 ] , op( map( op, [ seq( [ phi [ i ] , theta[ i ] ] , 
i = 2. .4 ) ] ) ) ] : 



# sum of mutual distances 

f := add( add( sqrt ( ( cos( Vt [ i, 2 ] )*cos( Vt [ i, 3 ] )- 
cos( Vt[ j, 2 ] )*cos( Vt[ j, 3 ] ) )~2+ 
( cos( Vt[ i, 2 ] )*sin( Vt [ i, 3 ] )- 
cos( Vt[ j, 2 ] )*sin( Vt [ j , 3 ] ) ) ~2+ 

( sin( Vt[ i, 2 ] )-sin( Vt [ j, 2 ] ) )~2 ), j = i+1..5 ), i = 1..4 ): 
f := spf ( simplifyC f ) ): 

Fls := [ seq( op( i , f ) , i = 1 . . 10 ) ] : 
p := [A, B, C, D, E ] : 

q := [ seq( seq( [i, j],j=i+1..5),i=1..4)]: 

# ten distances 

F := tableC [ seq( cat( p[ q[ i ] [ 1 ] ] , p[ q[ i ] [ 2 ] ] ) = Fls [ i ] , 
i = 1..10 ) ] ): 



# list of derivatives of f 

Df := [ seq( diff( f, indet [ i]),i = 1..7)]: 

# Hessian matrix 

H := MatrixC 7, 7, ( i, j ) -> ( diff( f, indet [ i ], indet [ j ] ) ) ): 

# coordinate corresponding the bipyramid distribution of 5 points 
Thetabp := [ phi [ 1 ] = -l/3*Pi, phi [ 2 ] = l/3*Pi, theta [ 2 ] = Pi, 
phi[ 3 ] = 0, theta[ 3 ] = -l/2*Pi, phi [ 4 ] = 0, theta[ 4 ] = l/2*Pi ]: 

# distance sum corresponding the bipyramid distribution of 5 points 

fmax := simplifyC subs( Thetabp, f ) ): fmaxrat := rfulb( fmax, 'r' , '1' ): 

py := ( 4+4*sqrt( 2 ) )*sqrt( l-tt"2 )+4*sqrt( 2 )*sqrt( 1+tt ): 
pyd := diff( py, tt ): 

# angle between the line connecting the spherical center and 

# the vertex of the pyramid bottom, and the pyramid bottom 
pyang := arcsin( solve( pyd, tt ) ) : 
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# coordinate corresponding the pyramid distribution of 5 points 

Thetap := [ phi [ 1 ] = -2*pyang, phi [ 2 ] = Pi/2-pyang, theta[ 2 ] = Pi, 
phi [ 3 ] = -arcsinC sin( pyang )*cos( pyang ) ), 
theta[ 3 ] = -arccot( sin( pyang )~2/cos( pyang ) ), 
phi [ 4 ] = -arcsinC sin( pyang )*cos( pyang ) ), 
theta[ 4 ] = arccot( sin( pyang )~2/cos( pyang ) ) ] : 

# distance sum corresponding the pyramid distribution of 5 points 
pyfmax := subs( Thetap, f ): 

# radius of the domain first excluded near coordinates corresponding 

# the bipyramid distribution 
xi := l/377*Pi: 

# radius of the domain first excluded near coordinates corresponding 

# the pyramid distribution 
Chi := l/791*Pi: 

# radius of the domain excluded near coordinates corresponding 

# the bipyramid distribution 
eta := xi: 

# radius of the domain excluded near coordinates corresponding 

# the pyramid distribution 
zeta := chi : 

bpyva := [ seq( indet [ i ] = Evalr( rhs( Thetabp[ i ] ) + [ -eta, eta ] 
i = 1. .7 ) ] : 

pyva := [ seq( indet [ i ] = Evalr( rhs( Thetap[ i ] ) + [ -zeta, zeta ] 
i = 1. .7 ) ] : 

# one domain need to be checked 

varsrangel := [ phi [ 1 ] = [ -2*arccos( ( fmax-2 )/9/2 ), ], 

phi[ 2 ] = [ 0, l/2*Pi ], theta[ 2 ] = [ 0, Pi ], 

phi[ 3 ] = [ -l/2*Pi, l/2*Pi ], theta[ 3 ] = [ -Pi, ], 

phi[ 4 ] = [ -l/2*Pi, l/2*Pi ], theta[ 4 ] = [ 0, Pi ] ] : 

varsrangel := map( x -> lhs( x ) = Evalr( rhs( x ) ), varsrangel ): 

# another domain need to be checked 

varsrange2 := [ phi [ 1 ] = [ -2*arccos( ( fmax-2 )/9/2 ), ], 

phi[ 2 ] = [ -l/2*Pi, ], theta[ 2 ] = [ 0, Pi ] , 

phi[ 3 ] = [ 0, l/2*Pi ], theta[ 3 ] = [ -Pi, ], 

phi[ 4 ] = [ -l/2*Pi, ], theta[ 4 ] = [ 0, Pi ] ] : 

varsrange2 := map( x -> lhs( x ) = Evalr( rhs( x ) ), varsrange2 ): 

k := 3: 

# subdivide the first domain 
ckrangelsl := intsbdv( varsrangel, k ): 

# subdivide the second domain 
ckrangels2 := intsbdv( varsrange2, k ): 

# used to test whether or not 5 points are on the same half sphere 
signdis : = [] : 

for i from 1 to 4 do 

for j from 'if'( i = 1, 3, i+1 ) to 5 do 
temp : = [] : 

nv := normal vector ( [ Vtd[ i ], Vtd[ j ] ] ): 
for k from 1 to 5 do 

if k = i or k = j then next: fi: 

temp := [ op( temp ), add( nv [ m ]*Vtd[ k][m], m= 1..3 ) ]: 
od: 

signdis := [ op( signdis ), temp ] : 
od: 
od: 

# gradients at points 

gr := [ seq( grad( i), i=1..5)]: 

grcv : = [ -gr [ 1 , 2 ] , -gr [ 1 , 3 ] , -gr [ 2 , 2 ] , 

gr[ 2, 1 ]*Vtd[ 2, 3 ]-gr[ 2, 3 ]*Vtd[ 2, 1 ] , 

seq( gr[ i, 1 ]*Vtd[ i, 2 ]-gr[ i, 2 ]*Vtd[ i , 1 ] , i = 3 . . 5 ) , 
seq( gr[ i, 1 ]*Vtd[ i, 3 ]-gr[ i, 3 ]*Vtd[ i, 1 ], i = 3..5 ) ]: 

# a condition used to test whether there exists a maximal distribution 

# in a domain 

grcv := spf ( cmb( expand ( grcv ) ) ): 
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# method names used to test 

Methods := [ def ver ,def eig, nondef , gradient, mindis, secondlength, 
dis, derivative, secordder, pointc, halfsphere ]: 

# save values of these variables in the file "f ivepoints .data" 

# if it does not exist 

save indet, Vt, Vtd, f, F, Fls, Df , H, 

Thetabp, fmax, f maxrat , pyang, Thetap, pyfmax, 
xi, chi, eta, zeta, bpyva, pyva, 
varsrangel, varsrange2, ckrangelsl, ckrangels2, 
gr, grcv, signdis, Methods, "f ivepoints . data" : 

else 

# read values of values from the file "f ivepoints . data" 
read "f ivepoints . data" : 

fi: 

# generate the process of subdisivion 
subdivisionprocess : 

end: 

# simplify expressions so as to raise the efficiency of interval computation 
spf := proc( expr ) 

local x: 

if type( expr, '+' ) and hastype( expr, radical ) then 

'+'( op( procnameC [ op( expr ) ] ) ) ) 
elif type( expr, '+' ) then 

map( combine, 

factorC '+'( op( selectC x -> nops( x ) = 5, [ op( expr )])))))+ 

'+'( op( select ( X -> nops( x ) <> 5, [ op( expr ) ] ) ) ) 
elif type( expr, '*' ) and hastype( expr, radical ) then 

'*'( op( procnameC [ op( expr ) ] ) ) ) 
elif type( expr, Or( '*', trig, arctrig, constant ) ) then 

expr 

elif type ( expr , < - ' ) then 

'"'( procnameC opC 1, expr ) ), opC 2, expr ) ) 
elif typeC expr, list ) then 

mapC procname, expr ) 
else 

error "unrecognized expression type: 7,0", expr 

fi 
end: 

# combine terms in denominators 
cmb := procC expr ) 

local tm, tmp, x, i: 

if typeC expr, list ) then return mapC procname, expr ) 

elif typeC expr, Not( ) ) then return expr 

fi: 

tm : = [ op C expr ) ] : tmp : = : 
for i to nopsC tm ) do 

if tmp = then 

tmp : = [ tm [ i ] ] 

else 

tmp := selectremoveC x -> denomC x ) = denomC tm[ i ] ), tmp ): 
if tmpC 1 ] = [] then 

tmp : = [ op ( tmp [2] ), tm[i] ] 
else 

tmp : = [ op C tmp [ 2 ] ) , 

C numerC op( tmpC 1 ] ) ) + numerC tm[ i ] ) ) /denomC tm[i] ) ] 

fi 

fi 
od: 

convert C tmp , ' + ' ) 
end: 

# convert from one coordinate system to another 
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# s2d: from spherical system to Cartesiaii system 

# d2s : from Cartesian system to spherical system 

# Cartesian system: [ xl, x2, x3 ] 

# spherical system: [ r, phi, theta ] 

# ( r >= 0, -Pi/2 =< phi <= Pi/2, -Pi < theta <= Pi ) 
changecoord := proc( L::list, opt ) 

uses IntervalArithmetic : 

if nargs = 1 then return changecoord( L, s2d ): fi: 
if type( opt, identical ( s2d ) ) then 

return [ L[ 1 ]*cos( L[ 2 ] )*cos( L[ 3 ] ), 
L[ 1 ]*cos( L[ 2 ] )*sin( L[ 3 ] ), L[ 1 ]*sin( L[ 2 ] ) ]: 
elif type( opt, identical ( d2s ) ) then 

ifL[l] =0andL[2] =0andL[3] >0 then 

return [ L[ 3 ] , Pi/2, ] : 
elif L[l] = 0andL[2] = 0andL[3] <0 then 

return [ -L [ 3 ] , -Pi/2, ]: 
elif L[l] =0andL[2] =0 then 

return [ , , ] : 
elif L[l] >0andL[2] = 0andL[3] = then 

return [ L [ 1 ] , , ] : 
elif L[l] <0andL[2] = 0andL[3] = then 

return [ L [ 1 ] , , Pi ] : 
elif L[ 1 ] > then 

return [ sqrt( L[ 1 ]'2+L[ 2 ]-2+L[ 3 ] ~2 ), 

arctanC L[ 3 ]/sqrt( L[ 1 ]~2+L[ 2 ] "2 ) ) , 

arctanC L [ 2 ] /L [ 1 ] ) ] : 
elif L[ 1 ] < and L[ 2 ] >= then 

return [ sqrt( L[ 1 ]"2+L[ 2 ]"2+L[ 3 ] "2 ), 

arctanC L[ 3 ]/sqrt( L[ 1 ]-2+L[ 2 ] "2 ) ) , 

arctanC L [ 2 ] /L [ 1 ] )+Pi ]: 
elif L[ 1 ] < and L[ 2 ] < then 

return [ sqrtC L[ 1 ]"2+L[ 2 ]"2+L[ 3 ] "2 ), 

arctanC L[ 3 ]/sqrtC L[ 1 ]-2+L[ 2 ] "2 ) ) , 

arctanC L[ 2 ]/L[ 1 ] )-Pi ] : 
fi: 
fi: 
end: 

# Euclidean distances 
eucdis := procO 

local i : 

if nargs = 3 then 

if typeC args [ 3 ], identical C sphere ) ) then 

procnameC changecoordC args [ 1 ] ), changecoordC args [ 2 ] ) ) 
elif typeC args [ 3 ], identical C descartes ) ) then 

procname C args [ 1 ] , args [ 2 ] ) 
else 

error "invalid argument: 7,1", args [ 3 ] 
f i: 

elif nargs = 2 and typeC args [ 1 ], list ) and typeC args [ 2 ], list ) then 
sqrtC sumC C args [ 1 ][ i ] -args [ 2 ][ i ] )"2, i = 1..3 ) ) 

elif nargs = 2 and typeC args [ 1 ], list ) and 

typeC args [ 2 ], identical C sphere ) ) then 
eucdis C changecoordC args [ 1 ] ) ) 

elif C nargs = 2 and typeC args [ 1 ] , list ) and 

typeC args [ 2 ], identical C descartes ) ) ) or nargs = 1 then 
sqrtC sumC C args [ 1 ][ i ] )~2, i = 1..3 ) ) 

else 

error "invalid argument: 7,0", args 
f i: 
end: 

# gradients 

grad := procC n:: integer ) 
local i: 
global Vtd: 

addC mapC '*', Vtd[ n ]-Vtd[ i ], 
1/simplifyC eucdisC Vtd[ n ], Vtd[ i ] ) ) ), 
i in [ $1..C n-1 ), $C n+1 )..5 ] ): 
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end: 



# normal vector of the plane determined by two points and 

# the spherical center 

normal vector := proc( L::listlist ) 

[ L[ 1, 2 ]*L[ 2, 3 ]-L[ 2, 2 ] *L [ 1, 3 ], 
L[ 1, 3 ]*L[ 2, 1 ]-L[ 2, 3 ]*L[ 1, 1 ], 
L[l, 1]*L[2, 2]-L[l, 2]*L[2, 1] ] 

end: 



# use various methods to check a rectangular domain 
ischecked := proc( invls ) 

local _method, i, j, intm, x, y, z, df , iv, tmp: 

global H, Fls , gr, f, f maxrat , Df , Digits, Vtd, grcv, signdis: 

uses IntervalArithmetic : 

_method : = : intm : = : 

if nargs > 1 then 

for i from 2 to nargs do 

if type( args [ i ], identical ( method ) =■[ 

# use vertex matrices to determine the negative definiteness of 

# the interval Hessian matrix 
identical ( defver ), 

# use eigenvalues to determine the negative definiteness of 

# the interval Hessian matrix 
identical ( defeig ), 

# use eigenvalues to determine the nonnegative definiteness of 

# the interval Hessian matrix 
identical ( nondef ), 

# use gradients to determine there exists no maximum in the domain 
identical ( gradient ), 

# check that the distance of some two points is less then 

# a certain value 
identical ( mindis ), 

# check that the distance of A and B is not the second longest 
identical ( secondlength ), 

# use interval computation to show directly that 

# the upper bound of f in the domain is less then fmax 
identical ( totaldis ), 

# use interval computation to show directly that 

# one of the derivatives doesnot change siges in the domain 
identical ( derivative ) , 

# analyze second order derivatives to determine that 

# one of the derivatives doesnot change siges in the domain 
identical ( secordder ), 

# check that C is below E, which contradicts assumptions 
identical ( pointc ), 

# check that all points are on the same half sphere 

# in which case f couldnot obtain its maximum 
identical ( half sphere ) 

> ) 
then 

_method : = rhs ( args [ i ] ) : 
elif type( args [ i ], 'Matrix' ( square ) ) then 

intm : = args [ i ] : 
else error "unrecognized method:'/,!", args [ i ]: 
fi: 
od: 
fi: 
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if _method = then _method := totaldis: fi: 
if memberC _method, {defver, defeig, nondef} ) then 
if intm = then 

intm := Matrix( 7, 7, ( i, j ) -> Evalr( H[ i, j ], invls ), 
shape = symmetric ) : 
fi: 

if _method = defver then 

isdef( intm, method = vertex, query = negdef ) 
elif _method = defeig then 

isdef( intm, method = eigenvalue, query = negdef ) 
else 

isdef( intm, query = nonnegsemidef ) 
fi: 

elif _method = gradient then 
for i to nops( grcv ) do 

X := Evalr( grcv[ i ], invls ): 
ifx[l] >0orx[2] <0 then 

return true : 
fi: 
od: 

return false : 
elif _method = mindis then 
for i to nops( Fls ) do 

if EvalrC Fls [ i ], invls )[ 2 ] < ( 2/15 ) then return true: fi: 
od: 

return false : 
elif _method = secondlength then 
0: 

EvalrC Fls[ 1 ], invls ): 
[] : 

for i from 2 to nops( Fls ) do 

z := [ op( z ), EvalrC Fls [ i ], invls ) ]: 
if z[ -1 ] [ 1 ] > y [ 2 ] then 
if X = then x : = 1 : 
else return true : 
fi: 
fi: 
od: 

if memberC false, map( tmp -> evalb( tmp[ 2]<y[l]),z))= false 
then 

return true : 
else 

return false: 
f i: 

elif _method = totaldis then 

return evalb ( Evalr ( f , invls ) [ 2 ] < f maxrat ) : 
elif _method = derivative then 

for i to 7 do 

df := EvalrC Df [ i ], invls ): 

if df[ 1 ] > or df[ 2 ] < then return true: fi: 
od: 

return false : 
elif _method = secordder then 
if intm = then 

intm := MatrixC 7, 7, C i, j ) -> EvalrC H[ i, j ], invls ), 
shape = symmetric ) : 
fi: 

df := □: 
for i to 7 do 
iv := []: 
for j to 7 do 

if intm[ i, j ] [ 1 ] > then 
iv : = [ op C iv ) , 

[ IhsC invls [ j ] ) = rhsC invls [ j ] )[ 1 ], 
IhsC invls [ j ] ) = rhsC invls [ j ] )[ 2 ] ] ]: 
elif intm[ i, j ] [ 2 ] < then 
iv := [ opC iv ) , 

[ IhsC invls [ j ] ) = rhsC invls [ j ] )[ 2 ], 
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lhs( invls[ j ] ) = rhs( invls [ j ] )[ 1 ] ] ]: 
fi: 

iv := [ op( iv ), [ invls [ j ], invls [ j ] ] ]: 
od: 

X : = Evalr ( Df [ i ] , map2 ( op , 1 , iv ) ) : 
y := Evalr ( Df [ i ], map2( op, 2, iv ) ): 
if type( X, constant ) then : 
else X : = X [ 1 ] : 
fi: 

if type( y, constant ) then : 

else y := y[ 2 ] : 

fi: 

df := [ op( df ), [ X, y ] ] : 

if df [ -1 ] [ 1 ] > or df [ -1 ] [ 2 ] < then return true: fi: 
od: 

return false: 
elif _method = pointc then 

if subs( invls, phi [ 2 ] )[ 2 ] - subs( invls, phi [ 4 ] )[ 1 ] < then 

return true : 
else 

return false: 
fi: 

elif _method = half sphere then 
for i to nops( signdis ) do 
X : = : tmp : = true : 
for j to 3 do 

y := Evalr ( signdis [ i, j ], invls ): 
if type( y, constant ) then 
if y > then y := 1: 
elif y < then y := -1: 
else y := 0: 
fi: 
else 

if y[ 1 ] > then y := 1: 
elif y[ 2 ] < then y := -1: 
else y := 0: 
fi: 
fi: 

if y = then tmp := false: break: 
elif y <> and x = then x := y: 
elif y <> and x <> y then tmp := false: break: 
fi: 
od: 

if tmp then return true: fi: 
od: 

return false: 
fi: 
end: 



# use specified methods to check a rectangular domain, 

# if not successful, the domain is subdivided, 

# and then check subdomains recursively, 

# until maximal width of subdomains is less than some value, 
spchecked := proc( va::list ) 

local X, y, methods, dim, subint, i, cur, bpyck, pyck, tmp: 

# "notchecked" and "checkprocess" should be set to [] before the procedure 

# is first called, since they store domains cannot be checked successfully, 

# and the process of checking. 

global bpyva, pyva, notchecked, checkprocess: 
uses IntervalArithmetic : 

methods := [ totaldis ]: 

# The case when 5 points form a distribution close to a bipyramid or 

# a pyramid should be checked separately, 
bpyck : = true : pyck : = true : 

if nargs > 1 then 

for i from 2 to nargs do 

if type( args [ i ], identical ( notcheckbipyramid ) ) then 
bpyck := false: 
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elif type( args [ i ], identicaK notcheckpyramid ) ) then 

pyck := false: 
elif type( args [ i ], list ) then 

methods : = args [ i ] : 
else 

error "invalid argument: 7,1", args [ i ]: 
fi: 
od: 
fi: 

if bpyck then 
for i to 7 do 

if not contain ( 
od: 

if i > 7 then 

checkprocess := 

return true : 
fi: 

for i to 7 do 
if rhs( va[ i ] 
or rhs ( va [ i ] 
break 

fi 
od: 

if i <= 7 then 

bpyck := false 
f i: 
fi: 

if pyck then 
for i to 7 do 

if not contain ( 
od: 

if i > 7 then 

checkprocess := 

return true : 
fi: 

for i to 7 do 
if rhs( va[ i ] 
or rhs( va[ i ] 
break 

fi 
od: 

if i <= 7 then 

pyck := false 
fi: 
fi: 

for i from 1 to nops( methods ) do 

# check the domain through specific in turn, 

# the order arranged for the methods may influence efficiencies 
if ischeckedC va, method = methods [ i ] ) then 

checkprocess := [ op( checkprocess ), [ i ] ]: 
return true : 
fi: 
od: 

# subdivide the domain over the widest interval 
dim := maxwidthdim( va ): 

# when the critical value is set to 1/1000, 

# domainl_1105_1101_1100_1099 cannot be check successfully, 
if intwidthC va[ dim ] ) < ( 1/10000 ) then 

notchecked := [ op( notchecked ), va ] : 
checkprocess := [ op( checkprocess ), [ -4 ] ] : 
return false : 
fi: 

# record the position of subdivision 
checkprocess := [ op( checkprocess ), dim ]: 
subint := intsbdv( va, 2, dim ): 



rhs( bpyva[ i ] ), rhs( va[ i ] ) ) then break fi 
[ op( checkprocess ) , [ -1 ] ] : 

) [ 2 ] < rhs( bpyva[ i ] ) [ 1 ] 
) [ 1 ] > rhs ( bpyva [ i ] ) [ 2 ] then 



rhs( pyva[ i ] ), rhs( va[ i ] ) ) then break fi 
[ op( checkprocess ) , [ ] ] : 

) [ 2 ] < rhs( pyva[ i ] ) [ 1 ] 
) [ 1 ] > rhs ( py va [ i ] ) [ 2 ] then 
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cur : = true : 

for i from 1 to nops( subint ) do 

# check subdomains divided 

if procnameC subint [ i ], methods, 

'if'( bpyck, NULL, notcheckbipyramid ), 

'if'( pyck, NULL, notcheckpyramid ) ) = false then 
cur := false: 

fi: 
od: 

checkprocess := [ op( checkprocess ), [ 'if'( cur, -2, -3 ) ] ]: 
return cur : 
end: 



# generate a list which indicates the process of checking 
subdivisionprocess := proc( ) 

local t, tl, X, y: 

global sp : 

t := [ $1. .2187 ] : 

tl := X -> op( map( y->y=t, x)): 
sp : = t : 

sp := subsopC tl( [ 62, 158, 239, 863, 1102, 1105, 1106, 2114, 2132] 
sp ) : 

sp := subsopC tl( [ [ 1105,1101 ], [ 1106,861 ], [ 1106,834 ], 
[ 1106,1099 ], [1106,1100 ] ] ), sp ) : 

sp := subsopC tl( [ [ 1105, 1101, 1100], [1106, 834, 725], 
[1106, 834, 726]] ) , sp ) : 

sp := subsopC tl( [ [ 1106,834,725,1752 ], [ 1106,834,726,1507 ], 
[ 1106, 834, 726, 1750 ] ] ), sp ) : 
end: 



# find the subdomain subdivided in the checking process 

# corresponding to an integer or a list integer in varsrangel 
subdivisiondomain := proc( dmcur::Qr( integer, list( integer ) ) ) 

local d, indm2, i, t, v: 

uses IntervalArithmetic : 

d := [ ] : 

indm2 := false: 

for i from 2 to nargs do 

if type( args [ i ], list ) then 
d : = args [ i ] 

elif type( args [ i ], identical (dm2) ) then 
indm2 : = true 

else 

error "invalid argument: 7,1", args [ i ] 

f i 
od: 

if d = [ ] then 

t := 'if'(indm2, varsraiige2, varsrangel ): 
d : = map ( rhs , t ) : 
fi: 

if type( dmcur, integer ) then 

V := Vector( 7, convert( dmcur-1, base, 3 ) ): 
[ seq( 

[ d[i] [l]+intwidth(d[i] )/3*v[i] , d[i] [1] +intwidth(d [i] ) /3* (v [i] +1) 
i=1..7 ) ] 

elif type( dmcur, list ) and nops( dmcur ) = 1 then 

procnameC op( dmcur ), args[2..-l] ) 
elif dmcur = [ ] then 

t := 'if'(indm2, varsrange2, varsrangel ): 

map ( rhs , t ) 
else 

procnameC dmcur [ -1 ], procnameC dmcur [1 . . -2] , args [2.. -1] ), 
args [2.. -1] ) 

fi 
end: 



# determine the subdomain in which a point or a domain is contained 
inwhichdomain := procC Is:: list ) 
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local b, indm2, d, i, v, m, t: 
global sp : 

uses IntervalArithmetic : 

if Digits < 100 then Digits := 100 fi: 

b := [ ] : indin2 := false: 

for i from 2 to nargs do 

if type( args [ i ], list ) then 

b : = args [ i ] 
elif type( args [ i ], identical (dm2) ) then 

indm2 := true 
else 

error "invalid argument: 7,1", args [ i ] 

fi 
od: 

# subdomain corresponding the list b 

d := subdivisiondomainC b, 'if'( indm2, 'dm2', NULL ) ): 

if type( Is, list( Or( rational, 'interval') ) ) then 
if memberC false, zip( contain, d. Is ) ) then 

print ( 'the point/domain input is not in any domain subdivided in the \ 
checking process' ): 
return FAIL 
fi: 

V := [seq( 

floor ( ( 'if'( type( ls[ i ], 'interval' ), 
ls[i][l], ls[i] )-d[i][l] )*3/intwidth( d[i]) ), 
i=1..7 ) ]: 

# when ls[i] = d[i,2] 

V := map( t->'if'(t=3,2,t),v): 
m := add( v[ i ]*3-( i-1 ), i = 1..7 )+l: 

# When the domain is in varsrange2, there is only once subdivision, 
if indm2 then 

#print( cat( 'the point/domain input is in:', 
#cat( domain2, _, op( map( t -> (t,_), b ) ), m ) ) ): 
return [ domain2 , op ( b ) , m ] 
fi: 

if op( [ op( b ) , m ] , sp ) = m then 

#print( cat( 'the point/domain input is in:', 

#cat( domainl , _, op( map( t-> (t,_),b) ),m) ) ): 

return [ domainl, op( b ), m ] 

else 

procnameC Is, [ op( b ), m ] ) 
fi: 

elif type( Is, list( constant ) ) then 

procnameC map( 'Evalr/shake ' , Is ) , b ) 
elif type( Is [ 1 ], '=' ) then 

procnameC map( rhs. Is ) , b ) 
else 

error "invalid argument: 7,0", args 
fi: 
end: 

# search the stored data to determine how a point can be checked 
howchecked := proc(ls : : list) 

local X, y, indm2, methods, t, fn, pd, currentdomain, i, cur, lev: 

uses IntervalArithmetic: 

if type (Is [1], ' = ') then 

return procnameC map( rhs. Is ) ) 
elif memberCtrue, map( t->not type( t, rational ), Is )) then 

return procnameC convert ( evalf C Is ), ft2rat ) ) 
fi: 

# first suppose that Is is in varsrangel 
indm2 := false: 

# input coordinates reprensent a distribution close to the bipyramid 
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if not member( false, zip( contain, bpyva. Is ) ) then 

print( cat ( 'via the "negative definite" interval matrix method, we had \ 
first exclude a rectangular region whose midpoint correspond to the \ 
bipyramid configuration, and radius is ', convert ( eta, name ), ', and this \ 
point just belongs to the region. ' ) ) : 

return 

# input coordinates reprensent a distribution close to the pyramid 
elif not memberC false, zip( contain, pyva. Is ) ) then 

print( cat ( 'via the "nonnegative definite" interval matrix method, we \ 
had first exclude a rectangular region whose midpoint correspond to the \ 
pyramid configuration, and radius is ', convert( zeta, name ), ', and this \ 
point just belongs to the region. ' ) ) : 

return 

# input coordinates are not in the variable ranges checked 
elif memberC false, zip( contain, varsrangel. Is ) ) and 
memberC false, zip( contain, varsrange2. Is ) ) then 

print ('the point is not in the checked range (-Pi/2<=phi<=Pi/2 and \ 
-Pi<=theta<=Pi) ' ) : 
return 

# Is is contained in varsrange2. 

elif not memberC false, zip( contain, varsrange2. Is ) ) then 

indm2 : = true 
f i: 

# methods used in the checking process 

methods := convert ( [pointc , halfsphere, mindis, secondlength, totaldis, 
gradient, derivative] , table) : 

methods [-1] := bipyramid: methods [0] := pyramid: 
methods [-4] := toosamlldomain: methods [-3] := fail: 
methods [-2] := success: 

# check which domain is Is contained in 

fn := inwhichdomainCls, ' if ' ( indin2, 'dm2', NULL ) ): 

currentdomain := subdivisiondomainC fn[2..-l], ' if ' ( indm2, 'dm2', NULL )): 

# input the path of the directory which stores checking data 
pd := readstatC "Enter the checking files directory:" ): 

if pd = NULL then 

currentdirC "D : /program/maple/paper/5 points on a sphere/data" ) 
else 

currentdir ( pd ) 
f i: 

# read a stored checking file 
read cat ( op ( map ( t -> (t , ' / ' ) , 

[ seq( cat( op( map( t -> (t,_), fn[l..i] )[l..-2] ) ), 

i=l. .nops(fn)-l)] ) ), cat( op( map( t -> (t,_), fn )[l..-2] ), '.txt' ) ): 

# print information of the file 

print( 'the point input is(expressed in float degree): ' ): 

print ( evalfC ls*180/Pi ) ): 

print ( cat ( 'it is in: ', \ 
cat( 'if'( indm2, domain2, domainl ), _, op( map( t -> (t,_), fn[2..-l] )\ 
[1..-2] ) ) ) ): 

print ( cat ( 'it was used: ', convert ( tm, name ), \ 
' seconds to check this domain' ) ) : 

cur : = true : lev : = : 

for i to nops( checkprocess ) do 

if cur and type( checkprocess [i] , integer ) then 
if Is [ checkprocess [ i ] ] <= 
( currentdomain [ checkprocess [ i ] ] [1] + 
currentdomain [ checkprocess [ i ] ] [2] )/2 then 

t := 1: 
else 

t := 2: cur := false: 
fi: 

currentdomain := intsbdv( currentdomain, 2, checkprocess [ i ] )[ t ]: 
elif cur then 

print ( 'it is finally contained in domain (expressed in radian): ' ): 
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print ( currentdomain ) : 

print ( 'that is (expressed in float degree): ' ): 

#print( 'it is finally contained in domain(expressed in float \ 

degree) : ' ) : 

print ( evalf ( currentdomain*180/Pi ) ): 

print ( cat( 'this domain is checked by method: ', \ 

methods [ op( checkprocess [ i ] ) ] ) ): 
break: 
else 

if type( checkprocess [ i ], integer ) then lev := lev+1 
elif checkprocess [ i ] in [ [-3], [-2] ] then lev := lev-1 
fi: 

if lev = then cur := true fi: 
fi: 
od: 

return 
end: 

end module : 



45 



